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Prcfifif 


The  focus  of  this  woik  is  timely  visualization  of  tiuee-dimensional  data  sets  doived 


from  computatitmal  fluid  dynamics.  The  rendering  algorithm  uses  an  octree  to  efEktiaitly 
traverse  die  data  set  in  rendering  an  iso«itfaoe  and  is  i^Jidicable  to  finite-differrace  as  well 
as  spectral  method  data  sets.  The  other  thrust  of  this  research  is  investigating  the 
advantages  (if  any)  of  rendering  a  qwctral  metiiod  solutitm  versus  a  finiteHlifference  one. 
As  a  result,  majority  of  the  results  is  centered  about  steady-state  solutions  to  a  model 
diffiision-ccMivectitm  problem  vriiidt  involves  a  boundary  layer.  In  addition,  limited  results 
for  a  driven  cavity  type  problem  is  provided. 
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Abstcifit 

The  timely  visualization  of  diree-dimatsioaal  data  sets  and  the  advantages  of  using 
a  spectral  mediod  sdutitm  versus  a  finite-dififeimice  method  solution  in  rendering 
isosuifaces  is  described.  The  Beam-Warming  numerical  algorithm,  which  uses  im{dicit- 
qipnndmate-factorization,  is  used  to  general  die  steady-state  solutions  for  a  model 
diiiusion-ctmvecticxi  problon.  The  Qidiyshev  collocation  operator  is  used  to  evaluate  die 
right-hand  side  of  die  Beam-Warming  algorithm  fw  the  qiectral  solution.  Gimparing  die 
model  problem  results  with  the  exact  solution,  the  spectral  series  solution  is  truncated  to  die 
same  degree  of  accuracy  as  the  firtite-difEerence  for  cmnparistm  of  rendering  times.  The 
rendering  algorithm  employs  octrees  to  efOciendy  traverse  the  data  set  to  fit  die  isosurfaces. 
The  actual  fitting  of  polygons  to  the  isosurface  uses  the  marching  cubes  table  lode  up 
algoridim.  Widi  the  spectral  series  solutitm,  interval  math  is  investigated  for  guaranteed 
detection  of  isosurfaces  during  die  initial  octree  travarsal(s). 
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RENDERING  OF  THREE-DIMENSIONAL  DATA  SETS 
DERIVED  FROM 

nNITE-DIFFERENCE  AND  SPECTRAL  METHODS 

L  Introdnction 

The  basic  problem  addressed  in  this  thesis  is  the  timely  visualization  of  the  results 
of  Con^nitational  Fhud  Dynamks  (CFD)  analyses.  Specifically,  this  research  investigates 
the  tendering  advantages  offered  by  the  use  of  qrectral  m^ods  (SMs)  over 
finite-difference  methods  (Fluffs)  when  sohiti(»s  are  tendoed  as  isosurfaces. 
Three-dimensional  (3D)  data  sets  ate  goierated  for  a  model  convection-diffusitNi  problem 
using  a  traditional  FDM  and  CSiebysliev  collocation  SM.  The  Beam-Warming  algorithm, 
an  implicit,  approximate-factorization  procedure,  is  used  to  generate  the  steady-state, 
finite-difference  solutions  and  is  modified  to  generate  the  q^ectral  solutions.  The 
Qiebyshev  collocation  opetatcn  is  used  to  evaluate  the  right-hand  side  of  the 
Beam-Warming  algorithm  to  provide  the  spectral  solution  at  steady-state.  Cmnparing  the 
model  problem  results  with  exact  solution,  the  spectral-series  solution  is  truncated  to  the 
same  degree  of  accuracy  as  the  finite-difference  for  compatistm  of  tendering  times.  The 
rendering  alg(»ithm  employs  a  linear  octree  to  efficiently  traverse  the  data  set  to  fit  the 
isosurfaces.  With  the  spectral-series  sohititm,  interval  mathematics  is  investigated  fm 
detecting  isosurfaces  during  die  initial  octree  traversal  Limited  implementation  of  the 
ccmventional  Beam-Warming  algorithm  to  a  driven  cavity  is  available.  Rendering  of  the 
FDM  driven  cavity  over  a  uniform  grid  is  provided. 
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Cunent  constraining  factors  in  rendering  3D  FDM  data  sets  are  die  size  of  the  dm 

set  and  the  spacing  of  the  grid  nodes.  For  a  3D.  a  100  x  100  x  100  grid  (1  million  nodes) 

is  not  uncommon.  Typically,  the  grid  nodes  ate  packed  closer  together  in  regions  of  lat]^ 

flow  gradients,  e.g.,  the  boundary  layer,  wakes,  etc.  Along  with  the  transformation  of  the 

idiysical  domain  to  a  uniform  rectaagular  ctxiqHitational  domain,  this  results  in  a 

non-uniform  distribution  of  nodes.  Such  grids  pose  a  serious  challenge  to  any  type  of 

volume  rendering  of  the  data  (Kerlick,  1990:2). 

To  be  useful,  volume  visualizadtHi  techniques  must  o^  understandable  data 
representations,  quick  data  manipulations,  and  reasonably  fast  renderi^.  Scientific 
users  ^ouM  be  t^le  to  change  parameters  and  see  the  reailtant  image  instantly. 

Few  present  day  systems  ate  capable  of  this  type  of  performance.  (Elvins, 
1992:194) 

However,  another  class  of  CFD  algcxithm  -  SMs  ~  may  be  able  to  bypass  these 
rendering  factors.  SMs  differ  from  FDMs,  since  they  provide  a  closed-form, 
truncated-series  solutitHi.  Typical  series  bases  are  the  Fourier  series  or  the  Chebyshev 
polynomials.  The  main  advantage  offered  by  a  SM  for  visualization  is  the  ability  to 
evaluate  the  closed-form  solution  where  needed  in  computing  an  isosutface.  The  series 
solution  can  also  be  highly  truncated  when  the  accuracy  level  required  for  volume 
visualization  is  much  lower  than  for  die  computation  of  the  solution.  Another  advantage  of 
SMs  is  the  abUity  to  use  interval  mathematics  in  evaluating  the  solution  over  a  grid  cell. 
Interval  methods  ate  guaranteed  not  to  miss  parts  of  contours  or  isosurfaces  down  to  a 
specified  size  in  the  viewing  region;  point  sampling  of  gridded  data  sets  cannot  guaramee 
this  (Suffem  and  FacketeU,  1991:331). 

In  addition  to  rendering  advantages,  SMs  in  general  require  fewer  degrees  of 
freedom  than  FDMs  (i.e.  less  nodes).  SMs  typically  offer  increased  accuracy  to  many 
otdersofmagnitucteover  a  FDM  for  the  same  nundier  of  grid  points.  Although  a  FDM  is 
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faster  than  a  spectral  method  whidi  uses  the  same  tamper  of  points,  dus  advantage  is 
negated  when  the  FEA4  uses  many  mme  grid  pdnts  to  produce  a  com|MBably  accurate 
sohidon  (Canuto,  et  aL.  1987:27).  The  final  advantage  is  the  compactness  (rf  data  sets  fOT 
comparably  accurate  sohitroos.  In  a  SM  series  solution,  any  individual  series  term  is  on  the 
order  of  the  term's  coefficient  Then  the  SM  data  set  can  be  fUtoed  to  drop  those 
coefficients  less  th«i  some  threshold,  i.e.  0.001 .  For  example,  when  the  model  problem 
SM  sohitioos  ate  truncated  to  the  same  order  of  accuracy  as  the  FDM  solutions,  the  ratio  of 
file  sizes  (SM:FDM)  is  roughly  5:7000  or  about  0.01  percent. 

Background  and  Prp^r  Work 

Presented  is  the  background  and  prior  wwk  on  the  numerical  and  tendering 
algorithms  used.  Since  the  intent  of  this  work  was  to  tender  a  3D  FDM  versus  a  3D  SM 
data  set,  the  Beam-Warming  numerical  algorithm  was  chosen  for  its  relative  ease  of 
implementatitHi  in  3D  and  ease  of  modification  to  a  Cfadtyshev  collocation  method.  FDM's 
basics,  the  Beam-Warming  algwithm,  and  extending  the  algorithm  to  SMs  ate  discussed. 
The  tendering  of  3D  data  s^  is  done  throu^  either  a  direct  vohime  tendering  or  a  throu^ 
asurface  fitting  algotithm.  A  surface  fitting  method  was  chosen  for  the  tendering  method 
since  these  m^hods  ate  faster  than  direct  volume  tendering  m^ods.  To  guarantee  surface 
detection  during  the  initial  traversal  of  the  data  set,  interval  mathonatics  is  used. 

12.1.  Finite-DifBeience  Basics 

Innite  difference  methods  provide  numerical  solutions  of  problems  in  fluid 
mechanks  based  upon  direct  apiaoximation  of  the  governing  equation(s).  For  example, 
the  sectmd-otder-accurale,  central-difference  oper^rs  derived  from  a  Taylor’s  series 
expansirm  ate  as  follows: 
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$  =  +  01  01  Ax^  h 

0  =  ****^:y.^4i  +OlAx^J  =  8^  +  OlAx^h 


(1.1) 

(1.2) 


where  ^xx  are  the  finite-difference  qteratois.  The  3  point  stencil  for  Eq  1  ^  is  shown 
in  Figure  1.1.  The  range  of  influence  is  local  and  encompasses  ±  1  additional  nodes.  As 
wiU  be  seen,  the  SM's  derivative  operators  have  global  influence  and  are  applied  over  tire 
whole  domain. 


Figure  1.1.  FDM  Three  Point  Stencil  in  ID 
1.2.2.  Beam-Wamring  Numerical  Algorithm 

Warming  and  Beam  (1978:85-129)  introduced  an  inq)]icit  apivoximate  factorization 
algorithm  for  a  system  of  conservation  laws,  such  as  the  cmiq^Tessible  Navier-Stokes 
equations.  The  algorithm  is  expressed  in  the  "delta  form"  (i.e.  solves  for  increments  of  the 
conserved  variable)  and  can  be  either  first- or  second-order  accurate  in  time.  Ifusedto 
march  to  steady  state,  the  accuracy  of  the  soludoo  is  governed  by  the  accuracy  of  the 
derivative  operators  and  boundary  conditions  applied  to  the  ri^t-hand  side  of  the 
algorithm.  The  algorithm  uses  the  Eulo*  and  viscous  flux  Jacobians  to  linearize  the 
Navier-Stokes  equations.  The  3D  flux  Jacobians  used  were  taken  from  Pulliam  and  Steger 
(1980:162).  By  spatially  factoring  the  finite-diiference  operator  on  the  left-hand  side  and 
by  neglecting  the  higher  order  terms,  the  discretization  of  the  left-side  in  3D  is  broken 
down  into  three  ID  system  of  equations.  By  doing  so,  and  since  the  ID  stencil 
encompasses  only  3  nodes,  the  algorithm  produces  a  block  tridiagonal  structure  which  can 
be  easily  inverted  at  each  time  step.  For  example,  Eqs  1.3  - 1.7  show  a  generic  algorithm 
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for  a  2D,  invisdd  flow  vriiich  is  factorized  into  two  1 D  systems  (die  Domeoclatme  a 
discussed  fiifther  in  Qiapters  2  and  3). 


[l  +  (1.3) 

[l  +5,A  +  5yB]il"(/  =  {[l  +M][i  ■^SyB\-S:,8yABY''V  (1.4; 

[i  +5,a+5,b]^"i;-[i  +m1[i 

[l  +M][i  +5yB}4*(/  =  9l"  (1.6) 

ri"I/  =  [l  +5yBp[l  (1.7) 


Without  approximate  factorization,  the  term  [l  +  *  bectxnes  conpaationally 

intractable  for  large  number  of  grid  nodes  and  fen*  the  extension  to  3D. 


1-2.3.  Sq££|q1] 


In  contrast,  SMs  differ  conceptually  frinn  FDMs.  Spectral  methods  belong  to  the 

class  of  weighted  residuals  methods.  Instead  of  directly  api»t>ximating  the  derivatives  as  is 

done  in  FDMs,  the  solution  is  first  assumed  to  be  the  form  of  a  truncated  series;  i.e.,  a  ID 

N 

Steady-State  solution  is  assumed  to  be:  Ifeie  ^j(x)'s  ate  known  as  the 

mO 

{q)proximating  or  "trail"  functions  and  a,'s  are  the  expansion  coefficients.  Then  die 
derivatives  are  tqiproximated  inditectly  by  die  derivatives  of  the  assumed  solution.  For 
example,  =  where  bj 's  can  be  determined  from  some  linear 

combination  of  a/s  and  from  the  orthogonal  properties  of  the  trail  functions,  s . 

Then  the  expansion  coefficients  are  determined  by  minimizing  the  wei^ted  residual  ~  the 
error  in  the  differential  equation  produced  by  using  the  truncated  approximation  instead  of 
the  exact  solution.  SMs  are  distinguished  by  the  choice  of  trail  functions  and  also  by  the 
type  of  weighting  method. 

The  choice  of  trail  functions  is  one  of  the  features  which  distinguishes  SMs  from 


other  weighted  residual  methods  and  FDMs  (Canuto,  et  al.,  1987:1).  In  SMs,  the  trail 
functions  are  defiled  ova-  the  whole  computational  domain,  are  infinitdiy  differmtiable. 


1-5 


and  are  orthogonal  to  each  other.  The  two  most  common  trail  fimctions  are  the 
trigonometric  and  the  Qiebyshev  polynomials.  For  evaluation  of  nonlinear  md 
non-constant  coefficient  terms  with  non-periodic  Ixmndary  condititms,  Gottlieb  and  Chszag 
(1977: 1 17)  recommend  using  Chebyshev  polynomials.  The  properties  of  the  Chebyshev 
polynomial  expansions  are  summarized  in  Appendix  C. 

In  determining  the  e3q>ansion  coefficients,  different  weighting  methods  can  be  used. 
The  3  most  ctnnmon  schemes  are  Galerldn,  collocation,  and  tau.  The  collocatitm  method  is 
the  simplest  to  implement  and  is  the  method  used  in  this  effort.  In  this  ai^roach,  the 
expansion  coefficients  are  determined  by  satisfying  the  differential  equation  exactly  at  the 
so-called  collocatitm  grid  nodes  (Canuto,  et  al.,  1987:1).  The  most  effective  choice  for  the 
grid  nodes  are  those  corresponding  to  quadrature  formulas  of  maximum  precision:  in  this 
effort,  the  Gauss-Lobatto  grid  nodes  were  used  (Canuto,  et  al.,  1987:13). 

Reddy  (1983)  discusses  the  application  of  pseudo-spectral  approximations  to  the 
evaluation  of  terms  in  the  Beam-Waiming  algwithm.  Pseudo-spectral  m^ods  ate  often 
interchanged  witii  collocation  methods.  However,  the  pseudo-spectral  evaluation  of 
nonlinear  terms  is  subject  to  ailiasing  errors  (Canuto,  et  al.,  1987:84-87).  Reddy  uses 
Fourier  pseudo-spectral  evaluation  of  the  right-hand  side  of  the  algorithm  and  also  uses  the 
thin  layer  Navier-Stokes  form  of  the  equations.  For  the  driven  cavity  jxoblem  solved  here, 
the  fill!  compressible  conservative  form  of  the  Navier-Stokes  equations  are  used. 

For  other  work  related  to  iq)proximate  factorization  schemes.  Street,  et  al., 
(1985:53-55)  uses  a  relaxation  scheme  in  which  tte  left-hand  side  terms  are  evaluated  with 
finite-difference  operator  and  the  right-hand  side  terms  are  evaluated  with  a  spectral 
'perator.  A  variation  of  their  scheme  is  used  in  this  effort. 
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I.2.4.  Renderiny  of  3D  nata  Sm 

The  nuy  (Mity  of  volume  vKualization  algoridims  are  designed  to  trader  gridded  3D 
data  sets  (i.e.,  scientific  and  medical).  These  volume  visualization  algorithms  fall  into  two 
categories:  direct  volume  rendering  and  surface  fitting  algorithms  (Elvins,  1992:196). 

1-2.5.  Direct  Volume  Rendering  Methods. 

These  methods  are  characterized  by  the  nuqtping  of  elements  direcdy  into  screen 
space  without  using  geometric  primitives  (e.g.,  triangles  or  quadrilaterals)  as  an 
intermediate  representation.  Instead,  each  volume  element  (voxel)  or  cell  contrUntfes  coltn 
to  the  final  image  through  a  data  classificatira  table.  Fdrexaiiq>le,  atid>lemaynuq>brae 
density  values  to  white/opaque,  muscle  density  values  to  red/semi-transpaient,  and  fat 
density  values  to  beige/mostly  tran^arent  colors.  Some  disadvantages  as  pointed  out  by 
Elvins  (1992:196-197)  are: 

•  to  get  a  reasonable  image,  care  must  be  taken  in  setting  up  the  classification 
table  ~  slight  changes  in  opacity  values  can  often  have  a  large,  unexpected 
impact  on  die  final  image; 

•  non-unifoim  grids  are  handled  with  difficulty. 

Specific  algorithms  include  ray  casting  (Foley,  et  al.,  1990:701),  Sabella  method  (Sabella, 
1988;  Foley,  et  al.,  1990:1037)  and  splatting  (Westover,  1990).  Since  direct  volume 
rendering  algorithms  are  more  computidonally  intrasive  than  surface  fitting  algorithms, 
surface  fitting  algorithms  are  considered  for  this  effort. 


1.2.6.  Surface  Fittiny 


Is. 


With  surface  fitting  algorithms,  the  3D  data  set  is  traversed  once  to  extract  a  surface 
of  constant  i^perty  value  or  an  "isosurface."  The  isosurface  is  fitted  with  geometric 
primitives  (e.g.,  triangles)  and  then  is  rendered  by  conventional  techniques.  Surface 
rendering  is  usually  faster  than  direct  volume  rendering  since  the  volume  is  traversed  only 
once  to  render  the  surface.  Some  disadvantages  as  pointed  out  by  Elvins  (1992:196)  are: 


•  changiiig  the  isosuiface  value  can  be  tune  consuming  since  the  endre  data  sM 
has  to  be  revisited  for  the  new  suifKe  extractkn; 

•  these  algorithms  suffer  problems  such  as  occasional  false  sinfaoe  pieces  and 
incorrect  handling  of  small  features  in  tire  data  set; 

•  lastly,  (mly  a  thin  surface  is  modeled  as  opposed  to  the  entire  data  set  in  direct 
volume  rendering. 

The  class  of  surface  fitting  algorithms  includes  contour  connecting,  opaque  cubes,  and 
matching  cubes.  In  addition,  an  octree  data-structure  can  be  inqdemented  to  reduce  the 
surface  extraction  and  re-extraction  times  (Wilhelms  and  Van  G^ler,  1992:205). 

Contour  conrrecting  was  one  of  first  methods  invented  for  volume  visualization 
(Elvins,  1992:197).  In  this  method,  closed  contours  are  traced  in  a  series  of  slices  and  then 
connected  to  form  an  isosurface  (Keppel  1975;  Fuchs,  et  al.,  1977).  The  two  main 
problems  widi  conrrecting  the  slices  are  die  "tiling  problem"  (how  to  best  graoate  the  mesh 
between  any  two  contours)  and  the  "branching  im^lem"  (how  to  best  tile  b^ween  slices 
which  contain  a  different  number  of  contours)  (Meyers,  et  al.,  1992:230).  Also  contours 
that  lie  within  another  contour  (i.e.,  a  cross  section  of  a  torus)  and  contours  from  two 
separate  surfaces  of  the  same  value  are  handled  with  difficulty. 

In  contrast  to  contour  connecting,  the  various  forms  of  the  cubes'  algorithm  take  a 
slightly  different  approach  in  fitting  the  surface.  Instead  of  dissecting  the  data  s^  into 
indivkhial  slices  and  thoi  fitting  the  geometric  jnimitives  to  the  slices,  these  algorithms  will 
fit  the  primitives  directly  to  the  data  set  cells.  The  oldest  cubes' algtnithm  is  the  qpaque 
binary  cube  or  cuberille  introduced  by  Herman  and  Liu  (1979).  This  algoridun  proceeds  in 
two  steps:  first,  the  volume  is  traversed  to  find  all  the  cells  which  straddle  the  isosurface 
value;  and  second,  the  cubes  are  rendered  to  ftnm  a  blod^  qrpearance  of  the  isosurface 
(Elvins,  1992:198). 

Lorensen  and  Oine  dranuttically  irrqnoves  this  widi  the  matching  cubes  algorithm. 
Instead  of  tendering  the  whole  cube,  one  to  four  triangles  per  cube  are  fitted  to  the  surface 
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inteisection  with  the  cube's  edges.  ThetriangulationisdtmediroaghatablekMdciqi 
covering  the  256  possible  cases  (Lotensen  and  dine,  1987:164-5).  These  256  cases  can 
be  bitdcen  down  into  14  distinct  nugor  cases  through  a  series  of  rotations  and  reflections. 
These  14  cases  are  shown  in  Hgure  1 2.  However,  of  these  14  mi^r  cases,  6  are 
ambiguous.  False  triangles  may  be  generated  if  special  handling  is  not  used  O^lhelms  and 
Van  Gelder,  1990;  Wyvill,  et  aL,  1986;  Nielson  and  Hanuum,  1991). 


Figure  1 2.  Matching  Cubes  Triangulation  Cases  (Locensen  and  Cline,  1987) 

As  a  means  to  reduce  the  surface  extraction  time,  Wilhelms  and  Van  Gelder  ( 1 992) 
introduce  branch-on-need  octrees  to  isolate  the  cells  straddling  the  isosutface  value.  In  a 
2D  planar  case,  a  quadrant  is  divided  into  4  squares  as  a  start  of  a  quadtree.  In  3D,  the 
volume  is  initially  divided  into  8  cubes  or  octants.  Then  each  octant  can  be  further 
subdivided  up  into  8  additional  octants  and  so  on.  Throu^  the  use  of  linear  octrees,  thra 
any  one  octant  can  be  miqiped  into  a  3D  array  containing  the  surface  values  (Gargantini, 
1981:367).  In  the  algorithm  used  by  Wilhelms  and  Van  Gelder,  an  octree  is  generated  with 
its  nodes  containing  minimum  and  iraximum  data  value  found  in  the  node's  subtree.  In  the 
surface  extraction  phase,  the  octree  is  traversed  with  a  particular  surface  threshold.  Only 
those  branches  containing  part  of  tire  isosurface  ate  visited.  If  an  octree  node  is 
encountered  with  its  maxunum  below  ot  its  minimum  above  the  threshold,  the  node  and  its 
children  nodes  are  ignored  since  the  octant  does  not  straddle  the  isosurface.  Whmi  a  node 
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at  the  bott(»n  of  the  tree  (a  kaf)  is  visited,  iwlygoos  are  generated  fipom  the  cells  makiiig  iq> 
the  leaf  by  using  the  marehing  ctdies’  table  look  iq>  (Wilhelms  and  Van  Gelda*,  1992:206). 
Since  isosuifaces  often  occupy  only  S  - 15%  of  die  total  vdume,  this  algorithm  provides  a 
means  of  quickly  sifting  through  the  volume.  They  report  surface  extraction  [diaae 
speedupsintherangeof  1.6to  11  and  ov«allspeedups(tf2-3over  die  regular  nuttchmg 
cubes. 

An  inqxxtant  time-saver  noted  by  Wilhelms  and  Vim  Gelder,  is  the  reuse  of 
conqxited  intmection  informadon.  Each  intersection  point  requites  interpolaiiiig  to 
generate  a  vertex  value  and  the  surface  normals  in  the  x,  y,  and  z  ditectioos.  They  save 
this  information  in  a  "hash  table.”  When  retrteved  for  the  last  dme,  the  infmnation  is 
deleted  to  make  room  for  new  vertex  infnmation.  The  branch-on-need  portion  of  their 
algorithm  efficiendy  allocates  memory  for  viewing  data  sets  whkdi  do  not  conform 
precisely  to  a  grid  dimension  of  2^x2^x  2^,  vriieie  ”d"  is  the  dimension  of  the  octree. 
However,  use  of  a  spectral-series  solution  negates  this  requirement  since  the  series  can  be 
evaluated  exacdy  at  the  octree's  nodes. 

Bloomendial  also  used  octrees  and  adiqidve  subdivisitni  to  organize  vertex  data  in 
Polygonization  of  implicit  functitxis.  However,  his  octrees  were  built  around  one 
isosurface  value  known  a  priori..  He  found  the  octree  in  this  use  a  convouent  mechanism 
for  adaptively  storing  the  implicit  surface  information  (Blomnenthal,  1988:354).  Also 
instead  of  using  Lorensen's  table-lookup  method  for  polygmiizing  the  cells,  he  used 
Wyvill's  (1986)  method.  The  advantage  Wilhelms  and  Van  Clelder's  algorithm  has  over 
this  is  the  abiliQr  to  change  the  isosurface  value  and  re-extract  the  resulting  surface  quicker. 

12...  Interval  Marhematica 

Suffem  and  Fackerell  (1991:331-340)  reports  on  using  interval  mathematics  in 
ctHijunction  with  octrees  in  rnidering  isosurfaces  of  implicit  functions.  In  the  above  octree 
algorithm,  a  surface  is  detected  if  any  two  nodes  of  a  cell  straddles  the  isosurface  value. 
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With  an  iiiq)licit  function,  die  cell's  range  of  vahies  can  be  used  to  calculate  the  natural 
iiuerval  extension  of  the  implicit  function.  This  natural  interval  extension  ainounts  to  die 
union  of  possible  variations  of  the  function  within  that  cell.  The  actual  functimi’s  value  is 
then  guaranteed  not  to  lie  outside  the  bounds  to  this  cell.  Thus,  inmval  mathemttics 
provides  a  100  percent  confidenoe  test  for  rejecting  octants  not  containing  the  desired 
isosurface.  The  main  disadvantage  of  interval  mathematics,  besides  requiring  an  implicit 
functi(Muil  descriptitm,  is  that  it  requires  additional  floating^ioint  operations  to  determine 
the  natural  extension.  Also  if  adiqitive  sub^vision  is  not  used,  then  the  usefulness  of 
interval  mathematics  is  constrained  by  die  final  {dot  depdi  or  the  final  resohitirai  of  die 
octree  algorithm  (the  surface  fitting  is  based  upon  the  cell's  actual  node  values,  not  its 
natural  interval  extension). 

1.3.  Summary  of  Resnits 

In  Chiqiter  2,  the  mathematical  formulation  of  the  model  problem  and  the  driven 
cavity  is  described;  Quqiter  3  •  the  numerical  algorithms;  duqiter  4  -  the  renderiitg 
algorithm;  Chapter  5  •  the  numerical  results;  and  Chapter  6  >  the  rendering  results. 

In  Cluqiter  3,  the  inqilementation  of  the  Beam-Warming  is  described  to  the  model 
problem  using  the  craiventional  second-oider-accurate  finite  difference  operators  and 
spectral  Chebyshev  collocation  operators  (m  the  right-hand  side  of  the  algorithm.  This 
inqilementatirai  is  for  grid  nodes  located  at  the  Gauss-Lobatto  grid  points. 

In  Quqiter  4,  die  modification  of  the  marching  cubes  algraidun  with  an  octree  is 
described.  The  application  of  interval  mathematics  in  conjunctitnwhh  a  ChdiyshevsNies 
solution  and  die  design  of  the  octree’s  attempted  hash  table  are  provided. 

In  Chiqiter  5,  a  compariscm  of  the  finite-diffetence  and  spectral  accuracy  to  the 
model  problem  is  provided.  Partial  impleiiientation  of  algorithm  to  the  driven  cavity  is 
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avail^le.  This  partial  implemeatttioii  is  for  finhe-difference  operators  on  a  uniformly 
space  grid.  Convergence  histories  of  the  algoriduns  are  i»ovided.  The  solutions  of  die 
numerical  algorithm  are  used  by  the  rendering  algorithm.  The  qiectral  data  sets  contaitiing 
the  expansion  coefficients  for  the  Ch^yshev  truncated  series  are  "filtered"  iniortouseby 
the  rendering  algorithm.  To  filter  the  model  problem  qiectral  data  s^  only  diose  disohite 
values  which  exceed  a  set  threshold  ate  retained.  The  resulting  accuracy  of  the  filter  data 
sets  are  on  the  order  of  10  times  the  threshold  used.  For  thresholds  of  0.001  -  0.00001, 
typically  80  -  250  expansion  coefficients  are  retained. 

In  Chapter  6,  tendering  times  fw  the  qiectial  and  finite-difference  model  proUem 
data  sets  are  cotrqiared.  The  initial  time  to  tender  the  first  isosutface  is  approximate  twice 
as  fast  for  a  qwctral  data  set  dum  for  a  finite-difference  data  set  of  conqiarable  accuracy. 
This  difference  in  timing  is  solely  a  function  of  the  relative  size  of  the  two  data  sets.  Even 
though  a  qiecttal  data  set  of  86  coefficients  tenders  twice  as  fast  as  finite-difference  of  33 ^ 
data  values,  a  finite-diffoenoe  data  set  of  65^  data  values  renders  slighdy  faster  a  spectral 
data  set  of  256  coefficients.  However,  (mce  the  initial  isosutface  is  tendered,  the 
regeneration  of  isosutfaces  for  the  two  methods  ate  equivalent.  In  successive  generation  of 
new  isosutfaces,  use  of  octrees  to  minimize  the  surface  re-extractimi  time  proved  to  be 
invaluable.  The  isosurface  build  time  with  octrees  are  found  to  be  2 -3  times  faster  than  a 
generic  marching  cubes  algorithm. 

In  addition  in  ChaptCT  6,  die  interval  mathematics  results  ate  provided.  In 
conjunction  with  a  Chdiyshev  series  solution  and  an  octree  traversal,  interval  mathematics 
were  found  to  have  limited  usefulness .  The  interval  mathematics  ate  used  in  calculating  an 
octant's  initial  minimum  and  maximim  values.  Except  for  the  leaves  located  at  the  bottom 
of  the  octree,  a  resulting  octant  range  of  vahies  is  found  to  be  too  large  for  a  series  solution. 
As  a  result,  the  entire  octree  is  essentially  traversed  to  tender  an  isosutface  -  at  which  point, 
the  generic  marching  cubes  algorithm  is  mcne  effective. 
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Also  in  Ouq^  6.  use  of  a  hash  tabk  in  OMijimctum  with  the  octiee  is  provided. 
The  hash  table  is  used  to  state  the  interpolated  isosuiface's  vertex  and  surface  ntninals 
values  for  neighboring  octant’s  use.  Using  the  design  discussed  in  Oiapter  4.  the  hash 
table  slowed  the  octree  traversal  by  a  factor  of  2  -  3,  thus  negating  the  octree's  ^peed 
advantage. 

Lastly  in  QuqMer  6,  the  quality  of  the  model  prc^lem  results  versos  the  exact 
soluti(»  and  some  driven  cavity  isosurfaces  are  provided.  The  model  problrai  isosurfaces 
d^rade  significantly  when  the  isosurface  value  (isovahie)  is  within  1  error  norm  of  the  data 
set’s  mariinmn  and  minimum  vahie.  The  maximum  error  norm  of  adata  set  is  cahailated 
by  OHnparison  with  the  exact  solution.  Driven  cavity  isosurfaces  are  provided  for  the 
pressure  field  and  die  velocity  magnitude  field  within  the  cavity 

Extrapolating  die  model  problem  results  to  die  tendering  of  a  qiectral  fluid  flow 
solution  is  not  straight  forward.  The  model  problem  frequency  content  inherently  is 
limited.  The  number  of  coefficients  required  for  a  fluid  flow  field  will  undoubtedly  be 
much  larger.  In  this  case,  the  FDM  will  have  lower  setup  costs.  The  one  advantage, 
which  has  not  been  fully  exploited,  is  the  ability  to  evaluate  the  SM  solution  at  any  arbitrary 
locatitm  in  die  donudn.  If  the  FDM  solution  is  defined  over  a  curvilinear  domain  (in  which 
the  node  spacing  is  a  function  of  two  or  rntne  coordinate  axes),  then  the  rendering  of  the 
FDM  solutions  would  require  interpolation  to  avoid  "cracks”  in  the  isosurfaces.  Then  the 
tendering  of  FDM  cases  will  significantly  inoease,  while  the  time  requited  fw  SM  cases 
would  remain  relatively  ctmstant.  This  area  warrants  further  investigatitm. 
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II.  M«th#madgal  ForialatfaMi 


Two  types  of  {MoUeins  aie  exanined  in  atempdog  to  obtain  q)ectnl  and  finite-differaioe 
datasets.  Tbefiist,  a  model  problem,  has  an  exact  solution  which  allows  comparison  of 
the  SM  and  FDM  absohne  accuracy.  The  secmnl,  die  driven  cavity,  was  chosen  as  a 
graetic  fluid  flow  and  for  ifsidative  ease  of  in^lmnemadon  over  the  dcmiain 
The  governing  equations  and  assumpdons  used  for  the  model  inoblmn  and  the  driven 
cavity  are  presented  herein.  Figure  2.1  shows  die  orientatkn  of  the  domain  and  die 
cavity’s  lid  boundary  condidon. 


Figure  2.1.  Model  Problem  ami  Driven  Cavity  Orientation 


2.1.  Model  Problem 

To  mimic  the  form  of  the  Navier-Stokes  equations,  the  following 
convectitm-diffusimi  equation  was  solved  numerically  for  the  model  problem: 

(2.1) 

where  Re-100.  Using  Diiichlet  boundary  conditions,  gfx,y,z)  and  boundary  data  are 
chosen  so  that  the  exact  sohititm  is 
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(2.2) 


f(x,  y.  z)  =  taab(X(l-  x^))8in(f  ( 1 + /;)siii(f  ( / + z^>) 

The  bypetbolk  tangent  functitm  was  chosen  in  order  to  miioic  aboundaiy-layer  in  the 
X  direction.  The  bouiulary-layer  thickness  is  then  controlled  by  A ,  as  listed  in  Table  2.1 
and  as  shown  in  Figure  2.2. 

Table  2.1.  Boundary-Layer  Thickness  for  Model  Problem 


A 

JPiL  1 1^ 

2 

0.75 

5 

0.31 

8 

0.18 

10 

0.14 

Figure  22.  Hypeibolic  Tangent  for  Model  Problem 


2.2.  Driven  Cavity 

For  the  driven  cavity  problem,  the  assumptions,  the  governing  equations,  and  the 
boundary  conditions  used  in  deriving  the  steady-state  solution  are  provided.  The 
steady-state  solution  is  drieimined  by  time-integrating  the  3D,  unsteady,  Navier-Stokes 
equations.  The  fluid  medium  is  assumed  to  be  compressible,  thereby  allowing  the 
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tiiiie-iiitegrati(m  to  be  cmied  out  with  die  efBcient  Beam-Wanniiig  algcmdun.  The 
conservative  form  of  the  equatimis  is  chosen  to  allow  for  inqilicit  treatmrat  of  the  inessuie 
term.  The  driven  cavity  is  defined  over  the  domain  Q,  as  discussed  tdiove  and  is  modeled 
subject  to  the  following  assumptions: 


•  thermally  and  calorically  perfect  gas  widi  the  specific  gas  constant,  y  =  14 

•  newtonian  fluid 

•  constant  specific  heats,  Cp  and  Cy 

•  the  lid  moves  with  acmistant  velocity,  U 

•  laminar  flow  with  Reynolds  number.  Re  =  100 

•  cmistant  walltenqieratute,  =  500K 

•  subsonic  flow  with  reference  Mach  number,  =0.3 

•  constant  Prandtl  number,  Pr  =  0.71 

•  no-slip  boundary  condition  at  wall 

•  nrmnal  derivative  of  pressure  at  the  wall,  ^  =  0,  where  r;  is  the  normal 
direction  to  the  wall 


The  reference  Mach  number  and  Reynokls  nurrdier  were  chosen  in  order  to  improve  the 
convergence  rate  of  the  Beam-Warming  algorithm.  The  constant  specific  heats,  gas 
constant,  and  Prandd  number  ate  standard  assun^itioiis  for  the  wall  temperature  given. 


22As  Naviflr^Srnkes  Rauafiomt  in  C.i 


itive  Form. 


The  governing  equations  ate  expressed  in  vector  form  as  (Pulliam  and  Steger, 
1980:159-160) 


31"^  dx'^  dy  dz  dx^  dy'^  dz’ 


(2.3) 


where  the  state  vector  of  conserved  variables,  V,  and  the  Euler  fluxes,  E,  F,  0,9X0 


defined  as 
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The  viscous  fhixes,  E,,  F,,  G^,  are  defined  as 
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(2.5a.b,c) 
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(2.6a,  b,c) 


Here,  p  is  the  density;  i^v,  ware  die  Caitesian  velocity  conqwnents;  p  is  die  pnessure;  T 


is  the  temperature;  is  the  total  energy;  e  istheintemalenagy;  k  is  the  thermal 
conductivity  coefBdoit;  are  the  shear  stress  componoits;  and  are  the  rangy 
components  related  to  die  shear  stresses  and  internal  raergy.  The  following  reladonships, 
derived  fimn  the  [neviousiy  mentioned  assun^itiras,  relate  the  diove  variables 
sqipropriately: 


p  =  (y-iM  (2.7a,b) 

E,^pe+jp(u^+v^+w^k  (2.8) 

T«  +V,  +Wj+2p^,  Ty  »  Tji  =p(^  +  ^),  (2.^b) 

A  +  +VT^  +  WTij.  (2.10) 
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2.2.2.  Boundary  Conditioiis 

On  the  statioiaiy  walls,  die  boundary  oHiditioiis  require  all  the  velocity  comiionents 
tovanish.  Density  is  related  to  total  energy  through  Eq  2.8.  The  final  boundary  andhions 
is  derived  from  the  assunqition  of  the  noniial  derivative  of  pressure  at  the  wall  to  be  zero. 

In  non-ditnensional  form  these  ate  as  follows: 
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fHimpVmpw^O, 

^  p » fijEf 

%-o  «  f-a 

where 

For  the  moving  lid,  Eqs  2.14a,  d,  e  are  modiiied  to  reflect  u  =  1: 

pu=p. 

pv-pwsO, 

whoe  ft  =[^ .[j 


(2.14a.b.c) 

(2.14d) 

(2.14e) 


(2.1Sa) 

(2.15b,c) 

(2.15d) 

(2.15e) 


IlL  Numerical  Alyorithm 


This  cluster  provides  an  overview  of  the  convmdonal  Beam-Warming  algorithm, 
including  an  analysis  of  a  variant  of  diis  alg(»itfam.  The  variation  incoiporates  die  eiqilidt 
Chebyshev  collocation  mediod  The  Chebyshev  collocation  procedure  provides  a  solution 
of  ^lectral  accuracy  at  steady-state.  The  algoridim  is  first  developed  widi  die  model 
problem  as  described  in  Chapter  3.  The  partial  extension  of  the  algorithm  as  described  here 
for  the  driven  cavity  is  implemented  for  finite-diffeience  on  an  unifmn  grid. 

The  accuracy  of  the  Beam-Warming  scheme  is  dependent  cm  die  choice  of  qiace 
and  time  differencing  operators.  Both  two-  and  three-time-level  sdiemes  may  be  used  with 
the  Beam-Warming  algorithm.  The  goieral  duee-level  scheme  can  be  written  as  follows 
(Warming  and  Beam,  1978:93): 

where  ri  is  the  fcnward-difiference  operates  and  V  is  the  backwaid-diffeieoce  operator 

A’^U  -  -  U*,  V"£/  =  !/•  - (3.2a, b) 

The  type  of  time  integratiem  represented  by  Eq  3.1  is  controlled  by  Oj  and  $2.  Table  3. 1 
lists  some  of  the  standard  schemes  and  their  accuracy. 


Table  3.1.  Typical  Integration  Schemes 


Integration  Scheme 

ej 

02 

Tisqiezoidal  Rule 

1/2 

0 

2nd 

Euler  Implicit 

1 

0 

1st 

3-Point  Formula 

1 

1/2 

2nd 

In  deriving  the  steady-state  solutions  within  this  effort,  die  Euler  implicit  integration 
scheme  was  used. 
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In  delta  form,  Eq  2.3  is  written  as 

-  A''E) + 4-(ArF,  -  A^P) + (3.3) 
at  ax  ay  Oz 

Through  the  introducticm  of  Hoc  flux  JacobUms,  Eq  3.1  may  be  written  solely  in  turns  of 

differences  of  C/,  A'^U.  The  flux  Jacolnans  are  used  to  linearize  die  flux  diffeienoes  in 

Eq  3.1  as  follows: 


where 


A’'EmA"A'*U  A’^F-^BTA^U  A’^G-C^A'^V,  (3.4a,b,c) 


(3.Sa,b,c) 


The  flux  Jacolnans  are  5  x  5  matrices  for  the  3D  driven  cavity  problem.  The  components 
of  all  the  Jacobians  matrices  are  given  in  Appendix  A.  The  viscous  terms  are  somewhat 
more  complicated  to  linearize  since  they  involve  duivadve  terms  of  the  state  vecttn  V .  The 
derivation  of  the  term  AflE^  is  shown  as  an  example  (Beran,  1993:AerD  753). 


(3.6) 

« jP"A"G  +  lf*A"l/,  =(F- +  (3-7) 

A"£;^  -(R"A''U)^.  (3.8) 

A’E,^  -Ar'E,^.  (3.9) 

A’E,-(R-A'‘V),+A’-‘E,^.  (3,10.) 

where 


The  term  (P-R^)  governs  the  spatial  variation  of  the  diermal  conductivity  and  viscosity. 


Since  this  term  will  ai^iear  cm  die  left-hand  side  of  die  algorithm,  it  has  no  effect  on  the 
steady-state  solution  of  the  driven  cavity.  By  assuming  die  viscosity  and  ctmductivity  are 
locally  constant,  the  turn  (F-  J^)  is  neglected.  Tte  term  AflE^^  pertains  to  the 
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cross-derivalive  variaUes.  It  is  handled  qybci^  in  order  to  pnxtoce  a  Mode  tridiagonal 
structure  for  the  algorithm.  The  odier  viscous  terms  are  as  follows: 


"-h 


where 


(3.12b.c) 


By  incorporating  Eqs  3.3  and  3.10  into  Eq  3.1  and  by  taking  advantage  of  apfwoximate 
factorization  yidds  the  Beam-Warming  algorithm  in  factored  form  (Beran,  1993:Aero  753). 

[/  +  ajfM"  “  +  atiSyB’'  -  5^")!/  +  )}d"£/  *  31".  (3.13) 

91-  =  a,[(£;  -  E)^  +(F,  -  FJ^  +(6^  -GJj,]"  (3.14) 

where  /  is  the  Identity  matrix,  and 

The  terms  Si,  Su  are  the  standard  FDM  derivative  operators  as  defined  in  Eqs  1.1  and  1.2. 
For  the  model  problem  described  in  Quqtto'  3,  Eqs  3.13  and  3.14  reduce  to  the  foUowing: 

[l*a,(S,r  + a,(S,f’  -lfc«=^]4V = 3f,  P'S) 

9t" =0.4""'/+ +/» +«, +«,j+«r  ■ 

where  e  =  and  g  is  defined  in  Eq  2.1. 

To  maintain  die  efficiency  of  the  algorithm,  the  boundary  ctmditions  must  be 
incorporated  while  preserving  die  tridiagonal  structure  and  achieving  the  desired 
steady-state  accuracy.  To  maintain  die  tridiagtmal  structure,  first-order-accurate  boundary 
conditions  are  implemented  for  the  implicit  portion  of  die  alg(»ithm.  Then  after  A’^U  is 
calculated  and  die  solution  has  bera  updated,  the  sectmd-tMtier  boundary  ermditions  are 
explicidy  enforced. 
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The  model  problem  is  implmnmiied  using  a  rectilinear  grid  comprised  (tf  inegulaiiy 
spaced  grid  points  at  die  Gauss-Lobatto  quadrature  points,  Xj  =  cos^.  The 
Gauss-Lobatto  points  are  the  extrema  of  the  Chebyshev  polynomials  (Gottlieb  and  Orszag, 
1977:40).  The  variation  of  the  grid  metrics  is  a  stretching/compressitm  of  any  one 
Cartesian  axis.  The  second-order-accurate  finite-difference  operators  are  modified  to 
account  for  the  non-uniform  grid  spacing,  Le.  =  x^+j  -  x,  (Anderson,  et  al.,  1984:55; 

Canuto,  et  al.,  1987:144).  For  the  boundary  nodes,  a  second-order  polynomial  fit  is  used. 
The  modified  operators  are  provided  in  ^ipendix  B. 

UL.  Beam  and  Warming  wi  Chebyshev  Collocation  Method 

At  steady  state,  the  solution's  accuracy  is  governed  by  the  method  of  evaluating  the 
right-hand  side  and  the  method  of  enforcing  of  the  boundary  conditions.  These  areas  are 
modified  to  yield  a  solution  of  spectral  accuracy.  Also,  a  relaxation  scheme  is  added  to  the 
algorithm. 

Instead  of  using  finite-difference  operators  to  evaluate  the  derivatives  of  the 
right-hand  side,  Eqs  3. 14  and  3. 16,  Chebyshev  collocation  differentiation  is  used.  This 
can  be  accomplished  through  the  use  of  a  Fast  Fourier  Transform  (FFT)  or  through  matrix 
multiplication. 

If  the  discrete  Chebyshev  transforms  are  computed  by  an  FFT  algorithm  which 
takes  advantage  of  die  reality  and  the  parity  of  the  function  u(  0)  =  u(cos0),  the 
total  number  of  operations  required  to  differentiate  in  physical  space  is 
( 5log2  N  +  8  +  2q)N,  where  q  is  the  order  of  the  derivative.  ...  If  the  collocation 
derivative  is  computed  by  matrix  multiplication,  the  total  number  of  operations  is 
2N^.  (Canuto,  et  al.,  1987:69-70) 

Canuto  et  al.  also  present  relative  timings  of  the  two  Chebyshev  methods.  These  are 
summarized  in  Table  3.2 
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Table  3.2.  Timing  (in  msec)  for  C^byshev  Collocati<»  Derivatives 


N 

Cyber  855  (scalar ) 

Matrix  Chebyshev 

Multiply  Transform 

Cyber  205  (vector) 

Matrix  Chdiyshev 

^Mt^^  Transform 

8 

0.11 

0.18 

0.002 

16 

0.38 

0.39 

0.004 

32 

1.54 

0.91 

0.009 

64 

5.74 

1.92 

0.046 

0.022 

128 

24.02 

4.20 

0.178 

0.044 

256 

93.49 

9.40 

0.706 

0.102 

(Canute,  et  al..  1987:45) 


In  vector  mode,  the  two  methods  are  equival^t  for  a  mediate  number  of  nodes,  Nj. 
Equivalency  is  achieved  between  N=16  and  N=32.  Also,  when  applying  the  FFT  method, 
ailiasing  can  be  introduced  in  the  evaluation  of  the  nonlinear  terms  if  care  is  not 
taken(Canuto,  et  al.,  1987:84).  One  method  is  to  use  SO  percent  m(»e  nodes  during  the 
FFT  (the  3/2  rule).  Then  the  additional  coefficients  are  set  to  zero  when  evaluating  die 
nonlinear  terms  via  the  convolution  sum  (Canuto,  et  aL,  1987:82-85).  For  the  spectral 
accuracy  required  here,  a  moderate  number  of  nodes  along  a  coordinate  direction  will  be 
used  (i.e.  25  -  33).  Since  the  transform  and  matrix  methods  require  iqiproximately  the 
same  amount  of  time  for  this  node  size,  die  matrix  method  was  used  to  avoid  the  issue  of 
ailiasing.  The  actual  derivative  matrices  are  listed  in  Appendix  C. 

As  for  the  boundary  conditions,  the  eiqilicit  updating  of  the  Neumann  conditions 
needs  to  be  modified  for  the  Chebyshev  matrix  derivative.  This  can  be  shown  for  the  ID 
case,  in  which  ^  =  0  at  n-l,N.  The  explicit  form  of  ^  is 


ai2 

022 


•••  oj/f 
•••  “2^ 

”•  ^NN. 


(3.17) 


where  A  =  is  the  derivative  matrix  and  17  s  is  the  array  of  coUocatirai  points. 
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At  bcmiKfanes,  die  fS(dkfwiiig  oonditk^ 


*  0  »  ujan  +  U2a2]-h  • 

^  =s  0  =  UiUjfi  +  U2a2S+‘  ’  +  “W«AW- 

Since  »  -onn  and  Cffj  =  -am  (aee  Ai^iendixC),  then  uj  and  Uff  can  be  calculated: 


(3.18) 

(3.19) 


“i  =  i(aiisj  +  affjSff),  (3.20) 

(3.21) 

where 

Sj  =U2a2i+"-+Ui^_jaff_jj, 

Sff  =  U2^fi-¥- . 

bsiaf/i-ajj. 

The  rdaxadon  sdmne  used  here  is  a  variation  of  the  scheme  used  by  Street  et  aL 
(1983:52-53).  As  they  pointed  out: 

The  crucial  property  dutt  a  relaxation  scheme  should  possess ...  is  that  it  damp 
effectively  tK  hi^-ftequency  components  of  the  error.  It  need  not  be  especially 
effective  in  die  low-frequency  range,  so  long  as  it  does  not  amplify  any 
components.  (Street,  et  al.,  1985:52) 

Their  modified  form  of  the  tqiproximate-factorization  algorithm,  expressed  for  two  qMce 
dimension,  is 


[or  -  S^rnial  -  SyH]^U  =  ocM,  (3.22) 

where  H  s  is  the  flux  Jacobian  of  M\  or  =  ^;  and  ate  the 

sec(xid-order-accurate  finite-differenoe  operates.  The  right-hand  side  matrix,  M,  is 
calculated  with  a  (Thdiyshev  collocadrm  method.  Their  relaxation  scheme  was  to  vary  a  in 
die  range  of  [a^.  by  using  the  rule 

«*  =  O*  W  .  (3.23) 

where  it  is  die  time-level  index  and  K  denotes  the  number  of  different  a  values.  Note 
that  Street's  a  starts  at  and  ends  at  Ui.  Since  Eqs  3.13  -  3.16  are  obtained  through 
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multiplication  die  following  adaptatkm  of  Eqs  3.22  and  3.23  is  used  for  diis 


investigation: 

(3.24a) 

where 

3x  ■  [a*f  ♦ 

(3.24b) 

3,  ■  [a*/  +  aj(SyB^  - 

(3.24c) 

3^  ■  [a*/ + a,(S,C^  - 

(3.24d) 

(3.25) 

and  k~1.2,....K.  Here,  the  a  starts  at  a/  andendsat  (Xj^.  Typical  values  used  are 
[0.1, 0.8  >  1.0]  for  [ai,ai^]  and  K  =  500- 1000. 


3.3.  Program  Overview  and  Flowchart 

Figure  3.1  presents  die  gaienl  flow  of  die  prognun  BW3D.  The  block  solver 
(Beran,  1993:ABn)  7S3)  is  vectorized  over  eadi  sweq>  line.  The  implementadons  of  the 
Chdiyshev  algoridims  are  udeen  from  Canuto  ^aL  (1987:^ppendix  B).  The  actual  q)ectral 
coefficients  are  extracted  from  die  converged  solution  using  Canuto's  ID  Qiebyshev 
transform  a{^died  in  3  successive  sweeps.  The  convergence  criteria  is: 

(3.26) 

where  £  is  set  to  (near  madiine  zero).  A  rdativdy  small  value  rtf  e  is  required  for 

the  spectral  solutions  to  converge  to  the  minimum  error  solution  as  determined  duou^ 
examinatitm  of  die  modd  problem  results. 

The  driven  cavity  outyut  for  graphics  consists  of  the  fdlowing  solution  param^ers: 
die  Carterian  vdodty  componmits,  u,  v,  w;  the  scalar  magnitude  of  die  velocity  vector 
field,  9  pressure,  p,  and  temperature,  r.  The  ovmall  objective  is  to 
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compare  the  rendering  (tfFDM  and  SM  solotkms  are  ctf  die  same  order  of  accuracy. 

To  achieve  this,  the  SM  solutions  are  "filtered*  by  using  (mly  those  expansion  coefiBcients 
whose  absolute  value  exceeds  a  set  threshold.  Typical  thresholds  are  0.001, 0.0001,  and 
0.00001.  By  doing  so,  the  size  of  the  SM  data  sets  were  typicaUy  reduced  to  die  mdo'  of 
100  coefficients.  The  resulting  accuracy  of  the  filtered  model  problem  data  s^  are 
tqiproximalely  10  times  the  direshold  used. 
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HguieS.l  Main  Flowchart  of  Numerical  Algoritfam 
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This  cha|>lBr  describes  hem  readering  of  a  singfe  param^er  data  set  is  accomplished 
widi  octrees  modified  mtalm  advantage  (tfiittsrvalmalfannatics.  The  im|dementtak» 
mechanics  for  multiple  parameters  (Le.,  the  driven  caviQr)  are  discussed  in  Chapter  6. 
Radker,  diis  chapter  discusses  the  implementation  of  the  generic  marching  cubes  algorithm, 
the  overiaying  of  the  octree,  the  traversal  algorithm,  die  attempted  hash  table  addition  (for 
reuse  of  vertex  data  among  octants),  and  data  structures  used  in  die  program. 


4^  Generic  Marching  Cnbes  ImplemciitntiQn 


The  marching  cubes  algmidim  can  be  divided  up  into  3  rendering  phases:  irqiut, 
extraction,  and  diqdav  T-^wt  covers  reading  the  data  sri  firom  disk  into  the  3D  data  array 
or  the  "Cube".  Exracdon  invtdves  traversing  every  cell  in  the  Cube  with  an  iscntafaoe 
value  and  fitung  triangles  to  those  cells  which  intersect  the  isosurfsoe.  Display  covers  the 
actual  the  process  of  sending  the  isosurfaoe  polygons  to  the  graiducs  hardware. 

4.1.1.  Input  Phase- 

For  large  data  sets  containing  samided  volume  data,  die  input  phase  takes  die 
majority  of  die  total  rendering  time  to  accomplish.  As  win  be  seen,  op  to  90%  of  the  total 
rime  can  be  qient  in  reading  an  ASCII  file  ftnmaL  In  the  case  of  filtered  SM  data  sris,  this 
time  is  minimal  since  the  file  siases  are  relatively  smaU.  hislead,  the  evaluadcHi  die 
truncated  spectral  series  at  every  node  locatkm  in  die  Cube  win  consume  die  majority  of  the 
time. 


£ 


Once  the  user  specifies  the  isosoiface  value  or  the  "dueshtdd",  every  cdl  of  die 
Cube  is  visited  by  swelling  in  pairs  of  jty  planes.  Ihe  algnidun  is  shown  in  Hgore  4.1. 


function  maicliing_cube8 


{ 


^  1st  xy  plane  of  Cube  data  and  ncxmals  at  the  nodes 

mteipolitojcy^dane 

for  (izsl;  iz<zdtnii;  iz-H-) 


{ 


} 


} 


pet  nextxy  plane  of  Cube  data  and  nonnals  at  the  nodes 

mterpolatejcy_jdane 

mlaiK>late_A.edgesJ)etweeQ_.xyj>lanes 

calculate  taUe  index  &  triangulate  polygons 

output  poison  data 

sw^xypli^ 


fiinctira  inleipolate_xy_plane 


{ 


for  (iy=0;  iy<ydim;  iy-H-) 
for  (ixs^,  ix<xdim;  ix-H-) 

{ 

if  (surface_cross_x jedge)  interpolate  x  values  and  nonnals 
if  (soiface_cross_yjedge)  intaix>late  y  values  and  normals 

1 


} 

define  surfaoe.crossjedge 

(node  1  value  ^  isovalue  and  iso'i^ue  <  node  2  value) 
(node  2  value  iso  value  and  iso  value  <  node  1  value) 


or 


Figure  4.1.  Matching  Cubes  Algmithm 


The  nonnals  at  the  Cube's  nodes  ate  calculated  by  normalizing  the  gradimit  of  parameter,  as 
shown  in  Eq  4.1.  The  gradient  is  approximated  ditough  a  first  order  FDM  along  each  of 
the  coordituue  axis. 

fx*  +  fyj + /z* 

-^JX  JZ 

(/i +/?+/!)* 

where  n  is  the  unit  normal  vector  of  die  paramo/.  If  a  cell  is  encountered  duu  intnsects 
the  isosurface,  polygons  are  generated  tepresmting  the  pordrni  of  isosutface  within  that 
celL  If  any  two  nodes  of  a  cell  straddle  die  isosurface  value,  dien  that  cell  intersects  die 
isosurface.  The  actual  ttianguladon  of  the  cell  is  done  throu^  Lmensen  and  Cline's  table 
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kxdoip  method  (1987:163).  For  cdls  containiiig  an  ambiguoiB  face,  the  emhiyrity  can  be 
lesfrived  by  evahiatmg  the  fimctum  at  the  ceatn'  (tf  the  ambiguous  face  CW^lhdms  and  Van 
Gelder,  1990:80).  However,  this  has  not  been  implemenied  fw  this  effort. 

4.1.3.  Display  Phase, 

In  this  effort,  instead  soiding  the  polygem  infonnatioo  directly  to  gr^>hics 
hardware,  the  information  is  stored  in  a  vertex  and  pdygon  array  for  later  diq>lay  at  the 
user  convenience.  This  allows  the  option  to  di^^y  multiple  isosurfaces  widiout  having  to 
traverse  the  octree  for  each  surface  legeitetatioiL  These  data  structures  are  briefly  described 
below  in  section  4.3. 

Octree  Implcmentotion 

Prior  to  discussing  the  ^recifics  of  die  rendoing  phases,  physical  layout  of  the 
octree  are  discussed.  The  rendering  phases  with  an  octree  are  similar  to  duu  of  marching 
cubes  except  for  setup  and  traversal  of  die  octree.  The  algorithm  will  be  described 
assuming  interval  mathematics  will  be  used  to  determine  the  initial  minimum  and  maximum 
values  of  the  octants.  Use  of  a  hash  table  is  attempted.  The  hash  table’s  design  is 
provided.  Lasdy,  incorporating  a  cutting  plane  into  the  algorithm  is  discussed. 

4.2.1.  The  Cube  and  die  Octree  Layout 

Hrst  the  Cube  (the  3D  data  array)  layout  is  described;  second,  the  octree  layout; 
lastly,  die  “dimension”  of  a  octree  leaf  is  related  to  die  number  of  cells  and  nodes  it  covers 
in  the  Cube. 

The  site  of  the  Cube  (die  number  of  nodes  ctmtained  within)  is  related  to  die  plot 
dqidi  resolution  by  the  formula:  ^2^  +  /]  x  ^2^  -f  x  ^2^  +  /j,  where  p  refers  to  the  plot 
dqith  of  the  algmithm.  For  this  effort,  the  maximum  plot  depth  used  is  6  (or  a  Cube  size 
of  65x63x65).  Higher  plot  depths  could  have  been  easily  used.  The  numba*  of  cells 
ctmtained  in  die  Cube  is  related  by  a  similar  formula:  ^2^]  x  ^2^]  x  ^2^]  =  8^.  For 
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example,  if  the  plot  depth  =  0,  that  the  Cube  would  contain  1  cell  and  8  nodes.  If  die  plot 
depth  =  1.  dien  the  Cube  would  contain  8  cells  and  27  nodes.  Also,  the  plot  depth  can  be 
thought  of  as  the  “octant  dimension”  of  the  Cube.  This  will  be  useful  in  relating  an  octree’s 
leaf  to  the  Cube’s  nodes  and  cells. 

Similaily,  the  dimension  of  the  octree,  odim,  is  defined  as  the  number  of  levels 
contained  in  the  octree.  The  top  level  of  die  octree  {pcdevel-  0)  coneqionds  to  die  entire 
Cube.  At  the  next  level  down  (ocdevel-  1),  the  Cube  is  subdivided  into  8  octants.  Then  at 
the  next  level  down  {ocdevel  s  2),  each  octant  is  subdivided  again  into  8  additional  octants, 
and  so  on.  The  number  of  octree  nodes  or  octants  at  any  one  level  is  detennined  by  the 
formula:  8*^^.  At  the  bottom  of  the  tree,  there  are  8**"  leaves. 

As  for  the  size  of  a  leaf,  the  number  of  Cube  nodes  and  cells  contained  within  a  leaf 
are  defined  by  similar  formulas:  #nodes  =  (2”  +  f)x(2"  +  f)x(2"  +  fj  and  #cells  =  8". 
Here,m  is  defined  as  the  octant  dimension  of  the  leaf .  For  sake  of  clarity,  the  dimension 
of  a  leaf  will  be  refened  to  as  the  malim.  As  noted  prior,  if  a  leaf  is  visited  which 
intersects  the  isosurface,  then  polygons  are  fitted  using  the  marching  cubes  algorithm, 
hence  the  mcdim  label.  Figure  4.2  shows  the  Cube  and  the  different  octree  levels  fiM- 
p  =  4,  odim  =  3,  and  mcdim  =  1.  In  this  example,  the  Cube  contains  17^  nodes  and  any 
one  leaf  ctmtains  8  cells  or  27  nodes  of  the  Cube.  The  plot  (^th  and  octree  dimension  and 
leaf  dimension  are  related  by:  p  =  odim  -i-  mcdim.  Also  note,  the  term  [2^  -i-  ij  is  the 
xdim,  ydim,  and  zdim  in  the  marching  cubes  algorithm.  Figure  4.1.  Following  Wilhelms 
and  Van  Gelder  suggestion,  the  minimum  leaf  dimension  was  mcdim  =  1  (1992:206).  If 
the  plot  depth,  p  =  mcdim,  and  odim  =  0,  then  the  generic  marching  cubes  algorithm  is 
implemented. 

Rendering  Phases  Overview. 

As  before  widi  die  marching  cubes  algorithm,  the  rmdeiing  can  be  divided  into 
input,  extraction  and  display  phases.  The  display  phase  is  the  same  as  befcxe.  The  amount 
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of  time  the  Other  two  phases  take  is  dependent  on  wh^her  interval  mathematics  is  used. 
The  impact  on  the  input  phase  will  be  discussed  next,  followed  by  a  discussion  of  interval 
mathematics,  surface  extraction  by  the  octree  tra^«rsal,  and  the  attem[Xed  use  of  hash  table 
for  temporary  storage  of  vertex  information  among  octants. 


Octlevel  -  0 


Ocdevel  - 1 


Octlevel  >  2 


Octlevel  >  3  Cube  Composite 

Figure  4.2.  Octree  Overlay,  Plot  Depth  =  4,  Octree  Dimensitm  =  3 


4.2.3.  Octree  Input  Kia.se. 

In  addition  to  allocating  and  assigning  data  values  to  the  Cube,  an  octree  is  created 
that  contains  at  each  node  the  maximum  and  minimum  data  values  found  in  that  node’s 
subtree  and  a  pointer  to  the  minimum  (x,y,z )  location  in  the  Cube.  Even  though  a  Cube's 
index  can  be  calculated  directly  from  die  octree  node’s  address  (Gargantini,  1981:366), 
Wilhelms  and  Van  Gelder’s  approach  is  adopted  to  speed  up  the  octree  traversal 
(1992:212).  Compared  to  the  generic  marching  cubes  input  [diase,  the  octree  input  phase  is 
slighdy  longer  since  the  octree  nodes  have  to  be  populated.  If  using  interval  mathematics  in 
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conjunctioa  with  a  SM  solutioo*  the  octree  nodes  are  assigned  minimim  and  maximum 
values  firom  the  tq>  down  tnnead  of  die  normal  bottmn’s  iq>  initialization.  Only  yihea  a 
leaf  is  visited  are  the  Cube  nodes  assigned  a  data  value. 

4.2.4.  Interval  Mediods. 

The  probability  of  detecting  contours  or  isosuifaoes  by  using  die  i»ocess  of  point 
sampling  (as  is  done  in  die  mandiing  cube  algorithm)  is  a  function  die  plot  dqith  used. 
Figure  4.3  shows  4  cases  in  which  point  sampling  will  fsdl  to  detect  a  contour  in  2D. 


(a)  (b)  (c)  (d) 

Hgure  4.3.  Undetected  Contour  Cases  (Suffem  and  Fackerell,  1991:332) 

Suffem  and  Fackerell  (1991:332)  also  provide  the  basic  operations  of  interval  mathmadcs 


which  are  as  follows: 

[a,b]+[c.d]^[a+c.b+dj,  (4.2a) 

[a.b]-[c.d]  =  [a-d,b-c],  (4.2b) 

[a,b]9  [c,dj  -  [min(ac.ad.bc,bd),max(ac,ad,bc,bd)J,  (4.2c) 

[a,b]/[c.d]^[a.b]9[j,jh  provided  0€[c,dj.  (4.2d) 


Eq  4.2d  can  also  be  extended  to  cases  when  0  €  /c,  J 7 ,  but  this  particular  operation  is  not 
tqiplicable  in  deriving  the  natural  interval  extension  for  the  tnincated  snies  used  here.  To 
derive  die  natural  interval  extension  of  a  function,  the  fundamental  property  of  intmval 
madiematics  is  used: 


ifxeX,  then/(x)€FfX),  (4.3) 

where  FfX)  is  calculated  by  substituting  the  interval  X=7a,h7  intoTfx)  and  using  die 
interval  operations  defined  in  Eqs  4.2a-d  (Suffon  and  Facknell,  1991:332).  In  applying 


,JI.J  w,ii. .  He J|.Mi -L  jJLpii, L,f^ 


diis  to  the  octree,  ev^  octant  is  defined  over  some  3D  inietvaL  By  calculating  the  qiectial 
series'  natural  interval  for  an  octant  yields  die  maximum  and  minimum  fm  diat  octant  The 
actual  surface  value  for  that  octant  is  guaranteed  to  lie  within  this  natural  interval  The 
natural  interval  extension  for  a  Chebyshev  polynominal  is  shown  in  Appendix  D.  This 
provides  a  means  of  branching  into  the  octree  without  first  evaluating  the  spectral  series  at 
every  node  location  in  die  Cube.  Since  the  surface  fitting  is  only  dependent  on  node  values 
contained  the  octree's  leaves,  the  interval  madiematics  is  not  applied  to  die  bottom  level  of 
the  octree.  The  primary  drawback  is  diat  every  natural  interval  calculatioi  for  an  octant 
requires  6  additional  series  evaluations  which  cannot  be  used  in  die  actual  surface  fitting.  If 
the  first  isosurface  selected  intmsects  a  small  portion  of  the  total  volume,  dioi  using  interval 
mathematics  is  a  good  trade-off.  Then  while  the  user  is  examining  the  first  few 
isosurfaces,  the  program  can  evaluate  die  remaining  unfilled  nodes  of  the  Cube  in  the 
background. 

4.2.5.  Octree  Traversal  and  Extraction  Pha.<e 
Every  node  in  the  octree  is  an  octant  which  is  made  up  of  eidier  8  additional  octants  or  (at 
the  bottom  level)  at  least  8  cells  or  at  least  27  data  values  in  the  CXibe.  Traversal  of  the 
octree  is  accomplished  by  perfimning  a  preordered  traversal  that  recursively  visits  each 
node  and  its  subtree.  Figure  4.4  shows  the  traversal  algorithm.  The  octree's  data  structure 
is  briefly  described  in  section  4.3.  The  one  modification  made  from  Wilhelms  and  Van 
Geldefs  setup  is  the  addition  of  octant  full  flag.  This  is  used  in  die  binary  decision  to  '*fill" 
the  octant  using  interval  mathematics.  The  “ijk”  index  is  the  precalculated  index  into  die 
Cube’s  array  and  stored  in  each  leaf. 
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functioii  Build^SurfiKe 

{ 

start  at  top  of  octree,  ocdevd  s  0.  get  Ae  top  octant 
while  (not  done) 

if(  octaitt  is  not  full )  fillOctant  using  Interval  Math 
if(  surfaceNoAcKnt )  get  Ae  next  octant 
ei« 

{ 

if(  octlevel  is  not  at  bottom  of  octree ) 

{ 

subdivide: 

drop  down  A  next  ocdevel  and  get  the  1st  of  8  octants 

} 

else 

{ 

Maidi  through  the  leaf  starting  at  ijk  mdex  of  Ae  ChAe 
get  the  next  octant 

} 

} 

while(  all  8  octants  at  this  level  have  been  visited ) 

{ 

if(  octlevel  is  back  at  Ae  top )  done = TRUE 
el» 

( 

go  back  vp  A  previous  octlevel 
get  Ae  next  octant 

} 

} 

} 

} 

define  surfaceNotPresent 

octant  min  >  isovalue  or  octant  max  <  isovalue 


Figure  4.4.  (Xxtee  Traversal  AlgoriAm  (Sufifem  and  Fackerell,  1991:339). 


4.2.6.  Ha.sh  Table  PeidgiL 

hi  addition,  a  hash  table  should  be  used  A  save  previously  calculated  surface 
vertices.  Since  each  surface  vertex  is  generally  used  in  four  neighboring  polygons,  having 
A  recalculaA  Ae  vertex  infcxmation  (the  location  an  edge  and  the  vertex's  ntxmal)  can 
negate  the  octree  speed  advant^e  (Wilhelms  and  Van  Gelder,  1992:221).  As  suggested  by 
Wilhelms  and  Van  Gelder,  each  cell  edge  in  the  Cube  is  given  unique  k^.  They  suggested 
using  Ae  leafs  index  inA  the  Cabc  as  the  basis  of  Ae  edge  key.  Specially,  Ais  is  given  as: 


4-8 


4  *  index  *  di9DCtion_code,  mime  directi(»_code  is  1  forx,  2  foty,  or  3  for  z  (Wilhdms 
and  Van  Gelder,  1992:223).  This  is  also  assuming  Ae  dimension  of  a  leaf  imabm)  is  1. 
Then  due  to  traversal  order,  the  3  edges  adjacent  to  the  ‘origin'*  of  the  cell  will  nevo-  be 
visited  again  (Wilhelms  and  Van  Gelda,  1992:213). 

In  the  design  attempted  here,  3  hash  ud)les  are  used,  one  for  eadi  coordinate 
direction.  This  allows  the  leaf  s  index  into  the  Cube  to  be  used  directly  as  the  key.  If  an 
edge  intersection  is  calculated  which  is  not  on  the  "(»igin  axis"  of  the  leaf,  dira  the  key,  die 
inteipolated  vertex  and  3  normals  values  are  sbxed  in  the  appropriate  hash  table  at  the  first 
available  table  entry.  If  a  leafs  edge  is  found  to  intersect  the  isosurface,  dien  the 
ai^ropriate  hash  table  is  checked  with  the  index  key.  If  the  edge  is  cm  a  "origin  axis"  of  the 
leaf,  a  delete  flag  is  set  for  that  hash  entry.  Widi  this  setup,  accessing  the  hash  tables, 
slowed  die  traversal  by  a  factor  of  2  -  3. 

4.2.7.  Modifying  the  ()ctiee  for  a  Cutting  Plane. 

If  displaying  multiple  isosurfaces,  a  desired  feature  is  a  cutting  plane.  A  simple 
cutting  plane  in  either  the  j[y,yz,  or  xz  plane  is  easily  implemented  with  an  octree.  Fot 
example,  if  yz  cutting  plane  is  to  be  set  at  x  =  0.3,  then  prior  to  marching  a  leaf,  the  x 
value  of  the  leaf  s  origin  is  checked.  If  it  is  less  than  0.3,  then  march  it  Otherwise, 
discard  it  and  move  on  to  die  next  octant.  The  modificadon  to  die  octree  traversal 
algorithm  is  shown  in  italicized  font  in  Figure  4.3. 
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if(  surfaoeNotPreseiit )  get  the  next  octtmt 

{ 

if(  octlevel  is  not  at  bottom  oi  octree ) 

{ 

subdivide: 

drop  down  to  next  octtevel  and  get  die  1st  (tf  8  octants 

else 

{ 

octant  X,  y,orz<  cutting  plane  x,  y,  or  z) 

March  the  octant  starting  at  ijk  locatUm  in  die  Cube 
get  die  next  octant 

) 


Figure  4.5.  Cutting  Hane  Algorithm  Modificatkm 


4.3.  Program  Data  Structures 

The  program  is  written  in  standard  C  programming  language.  The  data  structures 
for  4  of  the  key  variaUes,  die  Cube,  the  octree,  the  vertex  array,  and  the  polygcm  array,  are 
shown  in  Table  4.1. 


Table  4.1.  Key  Data  Structures 


Cube 

Octree 

Vertex 

Polygon 

intfiill; 

unsigned  intx; 
unsigned  inty; 
unsigned  int  z; 
float  nodeValue; 
vmtexn 
{  float  x; 
float  y; 
float  z; } 

int  full; 

unsigned  long  ijk; 
float  min; 
float  max; 

vCTtex  V 
{  float  x; 
float  y; 
float  z; } 
votexn 
{  float  x; 
float  y; 
float  z; } 

unsigned  long  vl; 
unsigned  long  v2; 
unsigned  long  v3; 
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V.  Nmierigai  lt« 


Presoited  heieiii  are  the  model  problem  accuracy  otw^wiaons  for  die  FDM  and  SM 
cases,  die  convergence  histmies,  die  filtering  of  the  SM  data  sets  and  their  resultmg 
accuracy.  In  this  context,  FDM  rqiresents  a  second-(»der-accurate,  finite-difference 
method.  SM  iqxesents  a  Qiebyshev  coDocadmi  method.  AH  results  pertain  to  steacfy-staie 
solutitms.  The  results  for  die  driven  cavity  at  die  Gauss-Lobatto  grid  nodes  are  not 
available.  Even  though  the  modification  of  the  finite-difference  derivatives  for  die  Gauss- 
Lobatto  nodes  wodeed  for  die  model  problem,  the  extension  to  the  driven  cavity  is  not 
working.  All  of  the  driven  cavity  results  are  firmn  using  FDMs  on  an  uniform  grid. 

Model  Problem  Accuracy  Comonrisons 

Each  numerical  solution  is  compared  with  the  exact  soludon,  and  the  absolute  error 
is  calculated  at  each  grid  node.  The  maximum  error  is  the  maximum  norm  of  the  absdute 
error  matrix,  i.e., 

^max  *  - ^adadaud  Lcr 

The  different  cases  analyzed  and  die  maximum  error  for  the  data  sets  are  shown  in 
Table  5.1  and  plotted  versus  the  number  of  degrees  of  fieedom  fot  a  ID  case  in 
Figure  5.1.  The  convogence  criteria  is  based  upon  the  maximum  norm  of  the  AV  being 
less  than  some  small  qisilon,  e  (see  Eq  3.26).  To  allow  for  the  SM  cases  to  converge  to 
their  lowest  error  solution,  £  is  set  to  near  machine  zero  (O  [10*^^]).  As  expected,  for  the 
same  number  of  degrees  of  fieedom,  SM  fvoduces  solutions  2-5  ordms  oi  magnitude 
more  accuralB  than  FDM. 
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The  overall  accuracy  of  die  SM  cases  is  a  function  oi  the  numbo’  of  nodes  inside 
the  boundary  layer.  The  number  of  Chebyshev  polynomials  required  resolve  boundvy 
layer  distances  is  on  die  order  where  ^  is  the  boundary  layer  thickness  (Gottlieb  and 
Orszag,  1977:40).  For  the  worse  case.  ^=0.14,  the  number  of  nodes  required  to  resolve 
tlus  is  3.  For  the  ID  case  in  which  the  degrees  of  freedom  is  17,  diere  are  just  3  nodes 
located  within  the  boundary  layer.  This  case  exhibits  die  worst  qiectral  error,  0.294.  In 
fact,  none  of  the  cases  in  which  Nsl7  would  converge  without  first  a  FDM  soludoiL  In  all 
odier  cases,  the  SM  would  converge  from  an  imiHilsive  start  ot  a  Unearly  interpcdmed  2D 
solution. 


Table  S.l.  Model  Problem  Arcuracy  for  FDM  and  SM 


A 

Boundary-Layer 
Thickness,  ^ 

N 

^max,  FDM 

W,SM 

2 

0.75 

17 

1.40x10*^ 

1.41x10"* 

25 

4.87x10*^ 

2.32x10*^ 

33 

2.59x10*^ 

2.27x10** 

5 

0.31 

17 

5.60x10** 

1.48x10*^ 

25 

7.28x10*^ 

3.11x10"* 

33 

3.73x10*^ 

8.90x10*** 

8 

0.18 

17 

no  ccmvergence 

4.50x10*^ 

25 

1.21x10** 

1.84x10*^ 

33 

5.71x10*^ 

1.80x10*^ 

49 

2.37x10*^ 

3.90x10*^ 

10 

0.14 

17 

no  ccmvergence 

2.94x10* 

25 

1.84x10** 

9.34x10*^ 

33 

7.70x10*^ 

3.89x10*^ 

49 

3.00x10*^ 

4.83x10*® 

65 

1.60x10*^ 

n/a 
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Figure  S.L  Solution  Accuracy  for  Numb^  of  ID  Degrees  of  Freedom 

2^  Converfcncc  History 

Typical  convergence  plots  are  shown  in  Figures  5.2  and  5.3  for  the  model  problem 
and  the  driven  cavity,  respectively.  The  model  prtrf)^  plots  are  firom  cases  in  uddch  a  2D 
solution  was  linearly  interpolated  to  the  3D  dranain  for  the  initial  start  The  SM  cases 
typically  required  a  much  smaller  initial  At  for  ctmvergence  Aan  die  FDM  cases  (At  ==  .01 
versus  .05).  For  the  model  problems,  a  simple  variable  time  step  is  implemented.  At  die 
end  of  each  iteration  step,  die  program  checks  the  rate  of  change  of  AI/  and  applies  a 
proportional  change  to  At  i.e.,  for  a  converging  solution.  At  is  increased.  The  change  in 
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At  is  also  weighted  to  favOT  a  decrease,  as  of^wsed  to  an  increase,  in  At  The  limits  of  At 
is  set  to  [0.001, 0.1]. 


In  the  cases  shown  in  Hgme  5.2,  the  initial  At  is  set  .01  and  .05.  With  die  initial 
start  conditions  described  {Hevious,  the  initial  AC/  ncmns  are  0.1  (v  less.  The  FDM  cases, 
in  general,  converge  relatively  quick  to  the  SM  cases.  After  a  cotain  norm  levd  is  reached 
tlKiugh  (in  diis  example,  die  FDM  cases  would  slow  dramatically  for  the 
machine-aox)  e  (Eq  3.26).  In  the  driven  cavity  cases  (Hgure  5.3),  convergence  is  defined 
as  e  =  0.001At .  Also  d»  initial  solution  is  an  impulsive  start  The  break  in  the  Ns:65 
case  is  a  result  of  performing  lOOiteratirmsat  AtsO.Ol  and  dien  restarting  the  code  with 
At  =  0.1.  If  the  initial  At  is  set  to  0.05,  die  code  does  not  converge  for  N  =  65. 

Filtering  of  Expansion  Coefficients  for  Graphics 

In  order  to  c<Mnpate  the  SM  and  FDM  solutions  at  the  same  d^ree  of  accuracy,  and 
in  order  to  reduce  the  computational  workload  for  evaluating  the  SM  smies  solutirm,  die 
SM  solutions  are  fUtered.  All  expansion  coefficients,  udiose  absolute  value  was  bdow  a 
certain  threshold,  are  discarded.  Various  thresholds  from  0.00001  to  0.001  are  used.  The 
effect  of  the  threshold  on  the  rendering  is  examined  in  (Thtqiter  6.  Figure  5.4  shows  die 
distributions  of  two  SM  data  sets  (A  s  10  and  N=25,33).  As  seen,  the  majority  of  the 
coefficients  are  within  the  interval  [-0.0001,0.(X)01].  Also  note,  as  the  number  of  degrees 
of  freedcMn  increases,  the  mode  of  distribution  shifts  toward  the  left  (10-^  to  lO-^O- 

The  resulting  accuracy  of  the  SM  data  sets,  aft»  filtering  at  0.001  and  0.0001,  are 
listed  in  Table  5.2.  If  the  initial  accuracy  of  the  data  set  is  less  dian  the  filtering  threshold, 
that  die  resulting  filtered  accuracy  is  ^roximately  1  (»der  of  magnitude  above  the 
threshold.  If  the  initial  accuracy  is  greater  dian  the  threshold,  die  resulting  filtoed  accuracy 
is  on  the  same  ordo'  as  the  initial  accuracy.  Also  ncxe,  the  number  of  coefficients  retained 
is  dramatically  reduced  from  the  origiiial  data  set  Since  the  data  sets  are  represmitatitnis  of 
an  analytic  function,  the  data  sets  should  be  inherendy  frequency  limited.  Hence,  the 
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reductioo  in  coe£Bcients  is  not  surprising.  Unfortunalely,  widiout  a  driven  cavity  sdution 
at  the  apprc^mate  grid  spacing  (ix.,  die  Gauss-Lobatto  poims),  the  extraction  of  die 
spectral  coefBcients  would  yield  impcactiadile  results. 


Hgure  S.4.  Histogram  of  Full  Expansion  Coefficient  Data  Sets:  >l  s  10,  N  =  25, 33 


Table  5.2.  Spectral  Accuracy  for  Hlieied  Data  Sets 


B 

N 

Full  Series 

^max 

Truncated  Series 

Threshold 

=  0.0010 

Threshold 

=  0.0001 

^max 

23 

%FUlr 

23 

%Filtr 

2 

17 

1.41x10'^ 

1.10x10*^ 

7.29x10*^ 

83 

D 

25 

2.32x10'^ 

1.10x10*^ 

H 

7.29x10*^ 

83 

33 

2.27x10** 

1.10x10*^ 

45 

0.1 

7.29x10*^ 

83 

5 

17 

1.48x10*^ 

1.02x10*^ 

66 

1.43x10*^ 

267 

5.4 

25 

3.11x10"* 

7.50x10*^ 

66 

1.17x10*^ 

112 

0.7 

33 

8.90x10*^ 

7.49x10*^ 

66 

0.2 

1.10x10*^ 

112 

0.3 

8 

17 

4.50x10'^ 

1,98x10*^ 

90 

1.8 

4.29x10*^ 

500 

10.2 

1.84x10'^ 

8.91x10*^ 

83 

0.5 

2.46x10*^ 

155 

1.0 

1.80x10"* 

8.75x10*^ 

85 

0.2 

1.17x10*^ 

144 

0.4 

49 

3.90x10*^ 

8.74x10*^ 

85 

0.1 

1.16x10*^ 

144 

0.1 

10 

17 

2.94x10*^ 

2.27x10*^ 

297 

2.79x10*^ 

1485 

^9 

25 

9.34x10*^ 

7.09x10*^ 

89 

6.48x10*^ 

237 

^9 

33 

3.89x10*^ 

9.38x10*^ 

86 

1.91x10*^ 

155 

^9 

49 

4.83x10*** 

9.38x10*^ 

86 

0.1 

1.73x10*^ 

156 

0.1  1 

Lastly,  the  number  of  coefficients  retained  is  also  related  to  tbe  initial  data  set's 
accuracy.  For  example,  for  the  case  in  which  the  threshold  =  0.0001  and  A  =  10,  the 
number  of  coefficients  retained  decreases  firom  1485  to  155  as  the  numbo:  of  degrees  of 
freedom  increases  from  17  to  33.  This  is  the  same  trend  as  seen  in  the  histograms  and  tlu 
shift  of  the  mode  towards  the  left.  It  relates  to  the  level  of  noise  or  error  contained  in  each 
coefficient  For  the  cases  in  which  N  =  33  or  49,  the  accuracy  of  the  full  data  set  is  the 
same  order  of  magnitude  as  die  threshold  or  lower  than  the  threshold.  For  this  reascm,  die 
filtered  data  sets  with  N  =  33  are  ccMnpared  widi  the  FDM  data  sets  with  N  =  33,  and  the 
filtered  data  sets  with  N  =  49  with  the  FDM  data  sets  with  N  =  69.  To  visualize  die  effect 
of  filtering,  a  ID  case  is  shown  in  Figure  5.5  for  the  hyperbolic  tangent  ftmction. 


5-7 


X 

Hgme  5.5.  Resulting  Tnincated  Fit  of  Hyperbolic  Tangent,  Threshold  =  0.001 


In  the  ID  filtered  case,  the  x  diiecticm  expansicm  coefficients  are  extracted  from  the 
3D  SM  data  sets  (N=17, 25, 33)  for  y,  z  =:  0.  The  extracted  coefficients  are  then  filtered  at 
a  threshold  =  0.001,  resulting  in  reducing  die  number  of  coefficients  to  14, 13,  and  1 1 
respectively.  The  plots  versus  the  exact  soludtm  are  shown  enlarged  in  Figure  S.S.  As 
seen,  the  solution's  error  exhibits  itself  as  an  oscillatory  behavior  about  the  exact  soludon. 
Since  the  Chebyshev  polynomials  are  closely  related  to  the  Fourier  cosine  series,  i.e., 
r„fcos0)  =  cos(n0),  this  behavior  is  consistait  with  a  truncated  Fourier  series  fit  to  a 
straight  line.  In  Chapter  6,  isosurfaces  are  shown  that  display  this  oscillatory  behavior  if 
the  desired  isosurface  value  is  within  1  error  norm  (e^)  of  the  data  s^’s  upper  limit 
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XL  RendwiM  Rgmite  —d  CobcImIom 


In  this  chi^jter,  the  following  results  ate  ixesented: 

•  rendering  timings  of  various  modd  {HoUem  cases, 

•  the  effect  of  die  data  set's  error  on  model  problem  isosurfaces,  and 

•  multiple  isosurface  views  of  both  the  model  problem  and  the  driven  cavity. 

The  dTect  of  interval  madiematics  and  die  attempted  use  of  a  hash  table  on  imidering  times 
is  discussed.  The  effect  of  using  a  cutting  plane  is  shown  in  die  multiple  isosutfaoes 
views.  The  cases  with  die  smallest  boundary  la^  are  used  primarily  for  the  comparistm 
discussion  since  these  cases  contained  die  largo*  maximum  error  norms,  (see  Eq  5. 1). 

The  model  problem  cases  used  are  listed  in  Table  6.1. 

Table  6.1.  Model  Problem  Data  Sets  Used 


File 

Original 
Deg  of 
Iieedom 

X 

Boundary 

Layer 

Filler 

Threshold 

#  Filtered 
Coef 

^mat 

FD3310 

33^ 

10 

0.14 

n/k 

n/a 

0.0770 

FD6510 

65^ 

10 

0.14 

n/a 

n/a 

0.0160 

SM3310 

33^ 

10 

0.14 

0.00100 

86 

0.0094 

SM4910.1 

49^ 

10 

0.14 

0.00010 

156 

0.0017 

SM4910.01 

49^ 

10 

0.14 

0.00001 

258 

0.0001 

SM338 

33^ 

8 

0.18 

0.00100 

85 

0.0088 

SM335 

33^ 

5 

0.31 

0.00100 

66 

0.0075 

SM332 

33^ 

2 

0.75 

0.00100 

45 

0.0110 

6.1.  Rendering  Timings 

The  various  tendoing  times  widiout  using  interval  mathematms  for  two  plot  depths  are 
listed  in  Table  6.2  (files  FD6S10,  SM3310,  SM4910.1  and  SM4910.01)  and  Table  6.3 
(files  FD3310  and  SM3310).  All  timings  are  from  a  Silicon  Graphics  VGX  R4(X)0, 
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SOMhz  workstation.  Columns  3  -  6  are  initialization  costs  in  setting  up  the  Cube  and  the 
octree.  Only  the  build  surface  timing  is  applicable  to  rendering  a  new  isosuiface.  With  an 
octree  at  the  highest  resolution  (plot  depth  =  6).  this  is  less  than  a  second.  Any  calculations 
of  normals  are  based  iqton  gradient  of  the  scalar  field  and  do  not  change  widi  subsequent 
new  renderings.  However,  the  decrease  in  calculating  the  normals  between  the  fnire 
marching  cubes  algorithm  (Octree  dimension,  odim  -  0)  and  die  octree  algoritiun  results 
from  calculating  the  norrruds  pertairting  to  that  isosutface  instead  of  the  whole  Cube. 


Table  6.2.  Rendering  Times  Without  Interval  Math,  Plot  Depth  =  6,  Isovalue  =  0.9S 


FUe 

Cube/Octree 
Alloc  &  Input 
(sec) 

Function  Eval 
&  Octree  Init 
(sec) 

Calc 

Nramals 

(sec) 

Bmld 

Surface 

(sec) 

0 

0.3 

2.6 

2.8 

0 

21.4 

2.6 

2.9 

0 

37.8 

2.6 

2.9 

0 

60.1 

2.6 

2.8 

1  Analytic 

0 

0.4 

3.4 

2.6 

2.9 

D 

51.9 

2.2 

54.1 

0.2 

MSM 

n 

0.4 

22.7 

23.1 

0.4 

B9 

0.5 

39.1 

39.6 

0.5 

WSM 

0.5 

62.2 

62.7 

0.4 

0.7 

Analytic 

n 

0.4 

4.8 

5.2 

0.4 

0.7 

FD6510 

5 

51.7 

3.4 

55.1 

MM 

0.4 

SM3310 

5 

24.4 

24.5 

0.7 

SM4910.1 

5 

40.7 

■11 

0.7 

SM49 10.01 

5 

0.5 

63.6 

64.2 

0.3 

0.7 

Analytic 

5 

0.5 

5.9 

6.4 

0.3 

0.7 

Interestingly,  the  build  surface  timings  for  the  FDM  are  approximately  half  that  of 
the  SM  data  sets  and  the  analytic  solutiotL  The  rraidering  algorithm  fits  approximately  half 
the  number  of  polygons  to  the  FDM  data  sets  as  compared  with  the  others  for  the  this 
isosurface.  For  instance,  the  number  of  polygons  for  the  FD6510  isosurface  case  in 
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Table  6.3.  Rendering  Hmes  ^tfaout  Interval  Math,  not  Dq)tfa  s  5,  bovahie  s  0.9S 


File 

Cube/Octree 
Alloc  &  Input 
(sec) 

Function  Eval 
ftOrtieelnit 
(sec) 

||mP|| 

BSI 

Calc 

Ntxmals 

(sec) 

Build 

Surface 

(sec) 

FD3310 

0 

6.4 

0.1 

6.5 

SM3310 

0 

0.1 

2.6 

2.7 

0.4 

Analytic 

0 

0.1 

0.4 

0.5 

0.3 

FD3310 

3 

0.2 

6.8 

0.1 

■9 

SM3310 

3 

2.8 

2.9 

0.1 

mm 

Analytic 

3 

0.1 

0.5 

0.6 

0.1 

■a 

FD3310 

D 

6.6 

0.3 

6.9 

0.0 

■SI 

SM3310 

U 

2.9 

3.1 

0.1 

mm 

■1 

0.1  1 

0.6 

0.7 

0.1 

Table  6.2  is  8303  versus  16295  and  16231  for  the  SM  and  analytic  solutions  re^)ectively. 
As  a  result,  the  time  to  build  the  isosurface  is  cut  in  half.  The  reduced  polygon  count  for 
the  FDM  is  a  result  of  the  node  spacing  used  in  the  Cube.  The  FDM  data  values  are  located 
at  the  Gauss-Lobatto  points  which  is  mote  sparsely  distributed  in  the  interior  of  die  Cube. 
As  a  result,  less  cells  are  traversed  in  the  surface  fitting.  In  contrast,  the  SM  and  analytic 
isosurfaces  are  extracted  from  a  Cube  and  octree  that  are  set  up  assuming  uniformly  placed 
nodes.  For  isosurface  values  which  are  located  in  the  denser  FDM  grid  region  (Le.,  an 
isovalue  =  0.1),  the  two  methods  produce  comparable  polygon  count  and  build  timings. 

In  comparing  the  octree  versus  the  pure  marching  cubes,  the  results  are  similar  to 
those  found  by  Wilhelms  and  Van  Gelder  (1992).  The  octree  significantly  reduces  the 
building  and  re-building  times  for  isosurfaces.  The  only  difference  hoe,  these  timings  are 
based  upon  not  using  a  hash  table  to  reuse  vertex  information  among  octants.  The 
attempted  addition  of  the  hash  table  is  discussed  below. 

In  comparing  SM  cases  with  FDM  ones,  die  SM  model  problem  cases  are 
approximate  twice  as  fast  as  die  FDM  cases  in  vdiich  the  accuracy  of  the  two  solutitxis  are 
comparable.  This  relative  timing,  for  the  domain  of  die  model  problem,  is  solely  a  function 


6-3 


of  either  reading  the  data  values  fixun  disk  (FDM)  or  evaluating  the  function  at  the  nodes. 
Since  the  model  problem  ftequency  ctxiient  is  inhermtly  limited,  the  number  of  coefficiems 
required  for  a  fluid  flow  fleld  will  undoubtedly  be  much  larger.  In  this  case,  the  FDM  will 
have  lower  setup  costs. 

The  one  advantage,  iK^h  has  not  been  fully  exploited,  is  the  ability  to  evaluate  the 
SM  solution  at  any  artntraty  location  in  the  domain.  If  the  FDM  solutimi  is  defined  over  a 
curvilinear  domain  (in  ii«^h  die  node  placing  is  a  function  of  two  or  more  coordinate 
axes),  then  the  rendering  of  the  FDM  solutions  would  require  interpolation  to  avend 
"cracks"  in  the  isosurfaces.  Then  the  rendering  of  FDM  cases  will  significandy  increase, 
while  the  time  required  for  SM  cases  would  remain  relatively  constant 

Another  attempt  to  eiqiloit  this  SM  advantage  is  in  the  use  of  interval  madiematics  in 
determining  an  octant's  initial  minimum  and  maximum  value.  The  reladve  timings  with 
interval  mathematics  are  shown  in  Table  6.4.  The "%  C^alculated"  column  is  the  percent  of 
the  C^ibe's  nodes  filled  during  the  octree  first  traversal.  The  Additional  Initializaticm" 
timings  result  from  filling  die  nodes  not  visited  during  that  first  traversal  and  can  be  dtme  in 
the  background  while  waiting  for  user  input 

Table  6.4.  Rendering  Times  With  Interval  Math,  Plot  Depth  =  6,  Isovalue  s  0.95 


File 

Odim 

w/o  Int  Math 

w/  Int  Math  | 

Total  Init 
BuUd 
(sec) 

Total  hut 
Build 
(sec) 

Addl 

hut 

(sec) 

Init 

%Calc 

SM3310 

4 

24.3 

27.2 

1.7 

100 

5 

25.9 

27.7 

7.4 

77 

SM4910.1 

4 

40.7 

42.8 

1.7 

100 

5 

41.6 

39.8 

10.8 

77 

SM4910.01 

4 

63.8 

1.5 

100 

5 

65.1 

16.1 

77 

SM332 

4 

13.4 

9.6 

63 

5 

14.5 

7.8 

gw 

39 
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In  this  i^lication.  interval  mathemalics  is  found  to  be  of  Imutod  use.  The  mam 
problem  is  the  namral  interval  extenskm  f(v  a  Chd)yshev  series  is  not  ^fective  until  the  4dt 
level  of  Ae  octree  (odims  S);  Ae  resulting  "width”  of  the  natural  interval  is  too  large  at  die 
higher  octree  levels.  The  widdi  of  an  interval  is  defined  by:  w([a,bj)mb-a.  If/fx)is 
continuous,  then  the  following  is  true  (Snyder,  1992:123): 

w(X)-^0  w(f(X))-^0  (6.1) 

For  instance,  the  natural  interval  extensitm  for  whole  domain  of  die  Cube  (the  tq;) 
octree  level)  for  SM3310  case  is  [-1.83,2.44].  The  widths  of  the  octant  interval  remain 
large  at  the  higher  octree  levels  to  include  all  isosuiface  values  and  die  whole  Chibe  is 
effectively  traversed  along  die  octree.  Octree  traversal  of  die  entire  Cube  is  relatively 
inefficient  as  compared  to  the  pure  marching  cubes  traversal.  Interval  madiematics  may  be 
useful  if  the  plot  depth  is  greater  than  6  or  if  recursive  subdivisimi  based  uptm  the 
isosurface  curvature  is  employed. 

Lastly,  die  hash  table  impact  on  build  times  is  discussed.  Table  6.5  lists  die 
increase  in  build  timings  for  whoi  the  hash  table,  as  described  in  Chiyiter  4,  is  enabled. 

The  number  of  hashes  refers  to  the  number  of  visits  made  to  the  table  and  isosurface  vertex 
information  was  found  from  a  previous  octant  As  shown,  using  the  hash  table  slowed  the 
build  tin  ^  by  a  factor  of  2  -  3.  Even  so,  the  octree  timings  without  the  hash  table  are 
significandy  faster  than  the  pure  marching  cubes  method. 

Table  6.5.  Hash  Table  Rendering  Times,  Isovalue  =  0.95 


File 

Plot 

Depth 

Odim 

Build 
w/o  Table 
(sec) 

Build 
w/ Table 
(sec) 

Increase 

Ratio 

# 

Hashes 

%  Vertex 
Reused 

Analytic 

6 

5 

0.69 

2.06 

3.0 

8522 

26 

SM3310 

6 

5 

0.65 

2.10 

3.2 

8614 

26 

Analytic 

5 

4 

0.17 

0.29 

1.7 

2151 

25 

SM3310 

5 

4 

0.17 

0.28 

1.6 

2152 

25 
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OnalitY  of  Numerical  Solution  Graphics 


Figures  6.1  -  6.S  show  isosurfaces  for  FD3310,  FD6S10,  SM3310,  SM4910.1, 
and  SM4910.01  data  sets.  For  these  isosurfaces,  the  accuracy  of  the  data  set  is  significant. 
i.e.,  the  isovalue  within  the  range  of  1  error  norm  of  the  data  set's  maximum  value.  For 
these  cases,  the  model  problem’s  maximum  value  is  1.  Each  figure  shows  the  degradation 
of  the  isosurfaces  as  the  isovalue  reproaches  1.  The  sequence  of  Hgures  6.1  -  6.5  shows 
the  subsequent  improvement  in  the  isosurfaces  as  die  accuracy  of  the  data  set  improves. 

In  all  cases,  whenever  die  isovalue  is  within  1  error  norm  of  the  upper  (or  lowm-) 
limit  of  the  set,  degradation  of  the  surface  occurs  (as  compared  to  the  tendering  of  the 
analytic  solution).  The  onset  of  die  degradation  of  the  FDM  isosurfaces  is  mote  gradual 
than  the  SM  cases.  For  instance,  FD3310  is  noticeably  degraded  for  an  isovalue  =  0.9638, 
while  1-  £max  -  0.923;  whereas,  SM3310  is  noticeably  degraded  for  an  isovalue  =  0.9923, 
while  1-  Emox  =  0.991.  Also  the  SM  isosurface  degradation  ^tead  relatively  uniform  over 
the  surface  and  is  oscillatory.  This  was  seen  in  the  ID  filtered  fits  to  the  hyi.erbo]ic  tangent 
function  in  Chiqiter  S.  Lastly,  Hgute  6.5  shows  by  decreasing  the  filter  threshold  to 
0.00001,  the  SM  data  set  essentially  tracks  the  analytic  solutitm  to  the  limits  of  the  display. 
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Figme  6.4.  Isosurfaoes  for  SM4910.1  Data  1-Ep^  =0.998 
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Hguie  6.5.  Isosuifaces  for  SM49 10.01  Data  Set  versus  Analytic  Solution 
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Hgure  6.6  pOTtnys  the  effect  of  viewing  multipfe  isosuiface  fw  the  model  proUem. 
A  cutting  plane  is  employed  to  eliminate  ponitms  of  the  outer  isosurfaoes  for  unobsmictive 
viewing  of  the  inner  isosurfaces.  Cunmitly,  these  surfaces  are  (^>aque.  However, 
ttanspaiency  could  be  easily  implemented  for  an  alternate  viewing  capalnlity. 

6.4.  Driven  Cavity 

Figures  6.7  and  6.8  show  various  isosurfaces  for  the  scalar  value  of  the  velocity 
vector  field.  Each  node  value  is  determine  by  the  fdlowing  relation: 
q  =  •(-  +  w^)^,  where  u,  v,  and  w  are  Cartesian  velocity  components.  The  vortex 

structure  is  easily  seeiL  In  both  figures,  die  xy  cutting  plane  has  been  used  at 
progressively  increasing  z  locations  for  isosurfaces  1  -2  or  1  >  5  reflectively.  In 
Figure  6.9,  the  isosurfaces  for  Ns65  case  has  relatively  smoodier  surface  than  for 
case.  This  "jagged"  effect  is  also  fiparent  in  the  isosurfaces  for  the  pressure  field,  as 
shown  in  Figure  6.10.  The  increased  accuracy  for  N=6S  is  readily  apparent  in  Figures  6.9 
and  6.10. 
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Figure  6.6.  Cutaway  Views  of  the  Model  Problem 
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Hguie  6.9.  Multiple  Pressure  Isosurfaces  for  the  Driven  Cavity 
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6.S.  Conclusions 


Extrqwlating  the  model  problem  results  to  the  rendering  of  a  ^)ectral  solution  of  a 
fluid  flow  is  not  straight  forward.  The  model  problon  frequency  content  is  inherently 
limited.  The  number  of  coefficients  required  for  a  flow  field  will  undoubtedly  be  much 
larger.  In  this  case,  the  FDM  will  have  lower  setup  costs.  The  one  advantage,  which  has 
not  been  fiiUy  exploited,  is  the  ability  to  evaluate  the  SM  solution  at  any  arbitrary  location  in 
the  domaiiL  If  the  FDM  solution  is  defined  over  a  curvilinear  domain  (in  which  the  node 
spacing  is  a  function  of  two  or  more  coordinate  axes),  then  the  rendering  of  the  FDM 
solutions  would  require  interpolation  to  avoid  "cracks”  in  the  isosurfaces.  Then  the 
rendering  of  FDM  cases  would  significantly  increase,  while  he  time  required  for  SM  cases 
would  remain  relatively  constant  This  area  warrants  further  investigation. 
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The  flux  Jacobians  aie  defiiied  as 


^dE 

_  dF 

B*-r- 

c-.#. 

(A.la,b,c) 

*dU 

dU 

dV 

dE,^ 

S«rr^ 

dGf 

(A.2a,b,c) 

3Uy 

dV, 

where  die  state  vecttv  of  conserved  variables,  I/,  die  Euler  fluxes,  £,  F,  G,  and  the 
components  of  the  viscous  fluxes,  ,  F,^ ,  ,  are  defined  in  Eqs  2.4  -  2.13.  Pulliam 

and  Steger  (1980:162)  provide  these  5xS  matrices  in  gen^  form  as  follows: 

0  Kt  Kc  0  ' 

d-Kj,(r-2)H  Ktu-Kj,(7-I)v  Kcu-Kj^(r~lh  h(r-l) 

Kb^^-vB  K^v-Kg(Y-l)u  B~Kg(Y-2)v  Kcv-Kg(r-l)w  K^ir-I) 

A.«,C=  Kj^w-Kc(y-1)u  Ktw-Kc(7-I)v  B-Kc(7-2)w  JEc(r-i)  ,  (^.3) 


where 

=  +  +  (A.4) 

e  =  K^u  +  KBV+KcW,  (A.5) 

For  a  uniform  grid  spacing,  to  obtain  A, 


Ka  =  1,  Kjf^Kc^O. 


(A.6) 


Inasunilarnuumer, 


R,S,T=^\ 

P 


-1^01+ ajKg) 
-v(a]-^asKf) 
-w(ai+asKG) 


[+as(u%+v%+w%] 


0 

0 

0 

Cj-t-ajKg 

0 

0 

0 

ai+asKg 

0 

0 

0 

aj+ajKo 

i^oj+ajKg)  v(a;+asKp)  w{a;+asKc)  H 


where 


«;=Tt  «5=7fc  «7=^. 

=M^+V^  +  W^, 

Likewise,  to  obtain  R  for  a  uniform  grid  spacing, 

Kg^l,  Kf^Kg  =  0. 


(A.7) 

(A.8) 

(A.9) 

(A.  10) 
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Appendix  B;  Fiaite-PifftrcBCt  EgmtigM  tai  GaBM -Lftliatte  Pointe 


For  the  Gauss-Lobatto  points,  x,-  =cos^,  i  =  0.1,...,N,^  second-order 

accurate  finite-difference  derivative  operators  are  as  follows: 

•  Central-Difference: 

2  2 

=  (A«<l«soii,ttal.l984:55).  (B  j, 

=  (C«i>uio,etA,1987:144),  (B-2) 

^*-7^ 


where 


fy-JtL-  h-r-T  h  -kthzl 


•  Forward-Difference: 


where 


K,:i+a) 


(R3) 


•  Backward-Difference: 


where 


S  f  -  +  fi-2 

hi.2(l+a) 


a  =  ^.  hi_i  =  x,_;-x, 
^-2 


(B.4) 


The  forward-  and  backward-difference  operators  are  derived  from  a  second-order 
polynomial  fit 
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AppaiUx  C!  Chebvihg¥  Pol¥ao«lai  Proyrti—  mmA  Parivativ  Matricaa 


The  following  properties  of  Qiebyshev  pdyncMnial  etpansions  are  taken  firom 
Gottlieb  and  Orszag,  (1977:Appendix  A): 


•  The  Qiebyshev  polynomial  of  di^ree  n,  T,(x),  is  defiaed  1^ 

TJcose)  =  cos(nd}.  (C.l) 


Thus,  Tq(x)  =  1,  r^fjcJssjc,  T2(x)~2x^ -1.  T2(x)-4x^  -3x,  and  so  on. 
•  The  Chd>yshev  polynomials  are  die  solutions  of  the  differential  equation 


/i  d  i,  2\V2dT. 


dx' 


(C.2) 


diat  are  bounded  at  x  =  ±7.  They  satisfy  the  orthoganality  relation 


lilT.(x)TJxl(l-x‘Y‘'‘dx  =  (C.3) 

where  Co~2,  c^~l  for  n>0,  and  s 0  for  n 9^ m,  iot  n-m. 

•  Some  properties  of  Chebyshev  polynomials  are 

Tn*l(x)  *  2xT^(x)  -  T^i(x),  (C.4) 

\Tn(x^^h  (C.5) 

T^(±1)  =  {±1)\  T2,(0)^(-ir,  T2n^j(0)^0.  (C6a,b.c) 

•  The  following  formulae  relate  the  expansimi  coefficients  a„  in  the  smies 

f(x)^iaj^(x)^  (C.7) 

H=0 

to  the  e^qiansion  coefficients  b„  of  the  series 

l/(x)=  ibJn(x),  (C.8) 

H=0 
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vdiere  c,  are  d^ned  in  Eq  C.3  and  fot  various  linear  (^xaton  L.  Two 
formulae  are: 


2  ^ 

d,*—  Ipo,, 

p+uoU 


ii.--  ip<p‘-H^)a,. 
dx‘  c,  ^+2  ' 

p^nvmk 


(C.9) 

(CIO) 


Canute  et  aL  (1987:69)  present  the  Chdiyshev  coUocatioii  derivative  operator,  Djy.  in 
matrix  fcmn  for  die  Gauss-Lobatto  points,  x,-  s  cos-^. 

J>Nf<Xt)=  i(Dg)^f(Xj)  (i  =  O.....N)  (C.I1) 

where 


(D«L= 


-X} 


(C.12a,b) 

(C.12c) 

-1 
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Awf adil  Dl  Natural  laterYai  E»tcario«  for  m  Chcbvhcv  FolTaoMdai 


In  order  to  determine  the  natural  interval  extrasum  for  a  Clid>ysliev  series,  die 
cyclic  behavior  or  the  ‘monotraki^  intervals*  of  the  ntfa<order  Chdiyshev  pdynomial  has 
to  be  taken  into  account 


For  example,  an  inclusion  function  evaluadcm  method  can  be  defined  for  the  cosine 
operator,  based  on  the  observation  that  the  cosine  function  is  monotonically 
decreasing  in  the  interval  [x2n,  x(2n+ 7)],  and  monotonically  increasing  in  the 
interval  [x(2n  +  l),K(2n -f  2j], for integm*  n.  (Snyder.  1992:123) 

Similariy,  a  method  for  evaluating  the  natural  interval  extmision  of  a  nth-order  Chdiyshev 

polynomial  IS  defined*  ^lecalhng  ^^*1,  ^  cos(n^)  x  ^  cos then 


T^(cose)^±l. 


(D.l) 


whennd  =  rrp,  p=:0,l,...,n.  Forexample,  the  x  positions  of  the  4th-orderChebyshev 


polynomial  extrema  are 

n^4,  d  =  -fp  (p^O, 12,3,4): 

0—0,  y,  x: 

x  =  l,  0.71,  0,  -0.71,  -1: 

T^  =  l,  -1,  1,  -1,  1. 


(D.2a-d) 


Givoi  the  following  interval  X-[a,b]  ^  0  =  [d^,,  0i,\  then  die  natural  interval 
extension  of  a  nth-order  Chebyshev  polynomial  is  defined  as 


Tn({ea.et\)* 

[-1,1\ 

-l,max  {T,(9,),T,(dt,)]\ 

.  [mn[T^{ej,T^(et,)},l\ 

\min{T^(ea).Tn(h)l 

[  max[T^(OJ,T^(et,)Y 


if 

if  andp^iseven,  p|,odd 

if  Pa  is  odd,  Pleven 

if  otherwise 


(D.3) 


where  p^  =  Pb  = 
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